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ABSTRACT 

The  recording  of  deep  crustal  reflections  at  near- vertical  inci¬ 
dence  along  a  90-kilometer  continuous  profile  in  southern  Alberta  has 
demonstrated  the  applicability  of  the  seismic  reflection  technique  for 
research  of  the  lower  crustal  section.  The  field  procedure  included 
the  use  of  multiple  holes  and  arrays  of  detectors  designed  to  discriminate 
against  the  shot -generated  long  period  surface  waves.  Velocity  or  fan- 
pass  filtering  is  shown  to  be  a  very  effective  means  of  enhancing 
individual  reflection  events. 

A  method  of  incorporating  the  effects  of  attenuation  into  syn¬ 
thetic  seismograms  has  been  developed.  Comparison  of  the  autopower 
spectra  of  selected  time  intervals  on  field  records  with  similar  spectra 
computed  from  theoretical  seismograms  i  s  an  effective  means  of  inverting 
the  seismic  wave  data  into  models  representing  the  properties  of  the  lower 
crust.  Such  studies  yield  information  on  the  variation  of  the  specific 
attenuation  factor,  Q,  with  depth  and  the  nature  of  the  transition  zones 
which  reflect  seismic  energy.  Thin  sills  of  alternating  high  and  low 
velocity  material  appear  to  satisfy  all  the  observational  data.  The 
layers  are  generally  less  than  200  meters  thick. 

From  an  interpretation  of  the  seismograms,  a  cross  section  of  the 
crust  in  southern  Alberta  is  derived.  Large  structural  variations  with 
relief  of  more  than  10  kilometers  over  horizontal  distances  of  50  kilo¬ 
meters  are  found.  On  the  basis  of  the  reflection  data  and  a  detailed 
refraction  survey,  good  evidence  is  given  for  major  steeply  dipping 
faults  in  the  middle  and  lower  crust.  The  interpretation  indicates  that 
the  M  discontinuity  has  similar  variations.  Combined  gravity,  magnetic. 
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refraction  and  reflection  data  provide  strong  evidence  for  a  buried  Pre- 
cambrian  rift.  From  the  gravity  and  magnetic  trends,  the  feature  has 
been  traced  for  several  hundred  kilometers  across  southern  Alberta  and 
under  the  Rocky  Mountains  into  British  Columbia. 
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CHAPTER  1 

INTRODUCTION 

1.1  Why  the  reflection  method? 

The  crust,  or  outermost  shell  of  the  earth,  constitutes  only  about 
1.5  percent  of  the  total  volume  and  one  percent  of  the  total  mass  of  the 
earth.  However,  this  relatively  thin  portion  of  the  planet  is  of  much 
greater  importance  to  our  understanding  of  the  earth’s  dynamics  and  his¬ 
torical  evolution  than  such  figures  might  suggest.  Beneath  the  earth's 
.crust,  whose  base  is  generally  defined  by  the  Mohorovicid  seismic  dis¬ 
continuity,  lies  the  upper  mantle.  The  crust  completely  encloses  this 
part  of  the  earth  and  all  regions  below  it — the  lower  mantle,  the  outer 
and  inner  cores.  Thus  the  direct  observation  of  these  deeper  domains  is 
prevented,  except  for  samples  which  are  derived  from  the  mantle  by  volcanic 
activity.  Consequently  we  must  rely  on  indirect  observations.  But  since 
these  measurements  are  necessarily  made  on  or  near  the  surface  of  the 
earth,  the  effects  of  the  crustal  section  are  included  in  the  observations. 
Before  unravelling  the  structure  of  the  deep  earth,  it  becomes  necessary 
to  know  to  what  extent  the  shallow  earth  has  affected  any  observations. 

In  many  cases  this  requires  detailed  knowledge  concerning  crustal  structure. 

As  evidenced  by  the  numerous  mountain  chains,  ocean  ridges  and 
deep-sea  trenches,  the  crust  itself  is  continually  undergoing  mammoth 
deformations.  Some  aspects  of  this  activity,  notably  earthquakes,  have 
proven  catastrophic.  By  studying  the  structure  of  the  earth's  crust,  can 
we  help  to  understand  the  underlying  causes  for  such  phenomena?  For 
example,  will  we  be  able  to  give  warnings  of  future  earthquakes?  The 
answer  is  hopefully  "yes".  But  this  requires  a  knowledge  of  not  just  the 
general  characteristics,  but  the  detailed  anatomy  of  the  crust,  as  well 
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as  the  upper  mantle. 

Of  more  economic  importance,  we  might  ask  what  role  the  tectonic 
evolution  of  the  crust  has  played  in  the  formation  and  distribution  of 
ore  and  mineral  deposits.  Are  there  any  correlations  between  deep-seated 
geological  phenomena  and  the  occurrence  of  known  mineral  reserves?  Once 
again  it  is  only  details  of  crustal  structure,  as  deduced  by  all  possible 
geophysical  methods,  which  can  help  to  answer  the  question. 

The  importance  of  the  aforementioned  and  other  considerations  have 
not  been  overlooked  by  the  scientific  community.  Such  cooperative  ven¬ 
tures  as  the  International  Geophysical  Years,  Project  VELA  UNIFORM  and 
the  Upper  Mantle  Project}have  provided  tremendous  impetus  to  geophysical 
investigations  of  the  earth's  crust  and  upper  mantle.  Much  has  been 
accomplished;  much  more  remains  to  be  done. 

The  geophysical  studies  of  crustal  structure  have  normally  involved 
one  or  more  of  the  following  methods  of  investigation:  refracted  body 
waves  and  surface  waves  from  earthquakes  and  from  controlled  explosions, 
gravity  and  magnetic  analyses,  electromagnetic  studies  both  of  natural 
and  of  induced  fields,  and  geochemical  observations.  Of  these  methods, 
refraction  seismology  using  controlled  explosions  has  proven  most  success¬ 
ful  in  delineating  gross  crustal  structure.  However,  the  wavelengths  of 
the  refracted  waves  limit  the  precision  of  the  method  for  detailing  the 
fine  structure.  As  well,  velocity  decreases  cannot  be  directly  observed, 
but  must  be  inferred  by  other  means.  Thus  none  of  the  foregoing  methods 
has  the  analytic  power  necessary  for  high  resolution  studies  of  the  deep 
crust . 

The  seismic  reflection  technique  is  singularly  suitable  for  such 
studies.  Observed  wavelengths  of  the  reflected  waves  are  several  times 
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shorter  than  their  refracted  counterparts,  thus  increasing  the  resolution 
of  the  method.  Common  exploration  seismology  techniques  such  as  continu¬ 
ous  depth  profiling  and  expanding  spreads  for  velocity  information  provide 
means  by  which  our  detailed  knowledge  of  the  deeper  crust  can  be  signifi¬ 
cantly  increased.  Spectral  analyses  of  waves  reflected  from  deep  horizons 
can  provide  information  regarding  the  nature  of  the  reflective  process. 

In  turn  this  may  increase  our  understanding  of  the  physical  processes 
which  can  occur  in  the  lower  crust.  By  combining  the  results  of  various 
analyses  with  those  from  theoretical  studies  it  should  prove  possible  to 
place  limits  on  elastic  properties  of  this  outer  shell.  Consequently,  the 
seismic  reflection  method  can  provide  a  powerful  tool  for  extending  our 
knowledge  of  crustal  structure. 

Some  of  the  expectations  for  the  reflection  technique  can  be 
attributed  to  advances  in  scientific  technology  over  the  past  decade. 
Magnetic  tape  recorders,  both  analog  and  digital,  provide  systems  with 
wideband  recording  and  a  dynamic  range  far  exceeding  that  of  the  photo¬ 
graphically  recording  apparatus.  The  data  may  then  be  played  back  through 
various  types  of  filters  to  enhance  the  desired  signals  relative  to  the 
unwanted  noise.  Digital  computers  have  greatly  facilitated  the  handling 
of  large  quantities  of  data  and  enable  the  implementation  of  sophisticated 
data  processing  techniques.  Together  these  facilities  make  the  effective 
separation  of  reflected  energy  on  the  seismogram  a  more  distinct  possibil¬ 
ity. 

But  the  reflection  method  is  not  without  its  limitations.  Con¬ 
siderable  effort  must  be  expended  in  eliminating  noise  and  surface  waves 
which  may  obscure  the  reflected  energy  arriving  at  near- vertical  incidence. 
It  is  also  necessary  to  cover  an  area  in  detail  if  miscorrelations  of 
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reflecting  horizons  are  to  be  avoided.  This  causes  the  method  to  be 
very  expensive  and  time  consuming.  However  by  intensively  studying  a 
few  small  areas,  then  combining  the  results  with  broad-scale  gravity, 
magnetic  and  refraction  surveys,  it  should  prove  possible  to  extrapolate 
the  results  over  a  widespread  area.  This  approach  has  been  followed  in 
the  present  research. 

1.2  History  of  the  reflection  method 

The  use  of  near- vertical  incidence  reflected  energy  as  a  means  of 
investigating  total  crustal  structure  is  a  quite  recent  innovation.  As 
far  as  the  writer  knows,  until  1968  the  University  of  Alberta  was  the 
only  North  American  institution  systematically  applying  deep  reflected 
waves  as  a  means  of  crustal  investigation.  The  situation  is  quite  dif¬ 
ferent  in  Europe  where  for  a  number  of  years  German  and  Russian  geophys¬ 
icists  have  been  using  the  method  to  advantage.  In  contrast  to  this 
fact,  for  more  than  35  years  seismic  reflection  prospecting  has  been 
successfully  used  in  North  America  to  discover  many  major  petroleum 
reserves  located  in  the  upper  sedimentary  strata.  Probably  the  first 
application  of  the  reflection  method  of  seismic  prospecting  took  place 
in  Oklahoma  on  June  4,  1921  (Schriever,  1952).  A  number  of  scientists — 
Dr.  Jo  C.  Karcher,  Dr.  W.  P.  Haseman,  Dr.  I.  Perrine  and  Mr.  W.  C,  Kite — 
from  the  University  of  Oklahoma  undertook  this  preliminary  experiment 
to  test  whether  the  fundamental  idea  had  any  validity.  Financial  prob¬ 
lems  slowed  the  introduction  of  the  method  to  the  general  industry,  but 
within  a  decade  it  had  achieved  recognition  as  a  viable  technique  for 
the  discovery  of  oil  and  gas  deposits.  No  one  involved  with  this  initial 
work  could  have  foreseen  the  success  which  was  to  follow  in  the  next  few 
decades.  Reflection  seismology  is  by  far  the  most  effective  and  widely 
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used  of  all  petroleum  prospecting  methods  (Tucker,  1968). 

While  development  of  the  reflection  technique  for  seismic  prospec¬ 
ting  has  continued,  the  extension  of  the  method  for  studying  crustal 
structure  has  been  slow.  The  first  published  report  of  reflections  from 
the  lower  crust  was  that  of  Junger  (1951).  By  extending  the  recording 
time  in  routine  seismic  prospecting  to  ten  seconds,  reflections  in  the 
time  range  of  7.0  to  8.5  seconds  were  recorded  in  Big  Horn  County, 
Montana.  On  the  basis  of  repeatable  results  and  identical  reflections 
for  one  shot  from  perpendicular  spreads,  Junger  effectively  argues  that 
the  energy  is  arriving  from  some  interface  in  the  deep  basement.  By 
assuming  a  velocity  distribution,  the  depth  to  this  horizon  is  placed 
in  the  range  from  18  to  21  kilometers.  Many  other  isolated  instances 
of  the  recording  of  deep  crustal  reflections  have  since  been  reported. 
Those  published  before  1961  are  well  documented  by  Steinhart  and  Meyer 
(1961).  In  a  similar  review,  James  and  Steinhart  (1966)  have  chronicled 
the  results  of  reflection  crustal  studies  from  1960  to  1965  in  many 
different  countries. 

Geophysicists  in  the  Soviet  Union  and  in  Germany  have  had  consid¬ 
erable  success  with  seismic  programs  which  include  the  recording  of  sub- 
critical  reflected  waves.  Before  considering  their  work,  a  few  more 
recent  cases  of  reported  deep  reflections  in  Europe  will  be  mentioned. 
Btth  and  Tryggvason  (1962)  reported  on  the  first  seismic  investigations 
of  crustal  structure  in  Fennoscandia.  Their  experiment  was  designed  to 
record  near- vertical  incidence  reflections  from  discontinuities  in  the 
crust  or  upper  mantle.  The  sources  of  energy  were  large  quarry  blasts 
(360  to  8000  kilograms  of  explosives)  and  the  receivers  were  individual 
vertical  geophones  whose  output  was  recorded  on  refraction  apparatus. 
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Considering  the  experience  of  our  own  work  and  that  of  other  investigators, 
it  is  not  surprising  that  with  their  method  the  results  were  inconclusive. 
Chamo  (1962)  reported  on  a  feasibility  study  to  determine  whether  the  re¬ 
flection  technique  could  be  used  to  investigate  deep  crustal  structure. 

The  results  were  positive  for  depths  up  to  15  to  20  kilometers  and  it  was 
suggested  that  further  work  on  the  procedure  in  the  Turkmen  S.S.R.  be 
carried  out. 

In  Greece,  Papazachos  et  al.  (1966)  have  reported  on  crustal  struc¬ 
ture  in  southeastern  Europe,  and  particularly  about  Athens.  Part  of  their 
experiment  included  recording  of  near-earthquakes  (A  =  0  to  80  kilometers) 
on  a  three- component  short-period  seismograph  station.  By  a  careful 
analysis  of  the  resulting  seismograms,  many  events  which  were  attributed 
to  reflections  from  within  the  crust  were  noted  and  the  squares  of  the 
travel  times  were  plotted  against  the  squares  of  the  distances.  Their 
resulting  curves  indicate  a  three- layered  crust,  an  interpretation  which 
agrees  with  their  conclusions  from  the  analysis  of  refracted  and  direct 
waves.  No  indication  is  given  as  to  how  many  events  may  have  been  picked 
that  wouldn't  fit  the  plots.  In  a  similar  study,  White  (1969)  reported 
on  some  unusual  seismic  phases  recorded  in  South  Australia.  From  a  small 
network  of  seismograph  stations  set  up  to  monitor  local  earthquakes,  a 
phase  reflected  from  the  Mohorovicid  discontinuity  was  regularly  recorded 
and  identifiable  out  to  a  distance  of  A  =  2°.  In  addition  some  converted 
phases  of  the  S  to  P  type  reflected  from  the  M  discontinuity  were  found. 

It  is  important  to  note  from  such  studies  that  with  good  instrumentation, 
reasonable  evidence  is  provided  for  the  recording  of  sub-critical  reflected 
arrivals  from  within  the  crust  when  the  source  is  a  shallow,  low-magnitude 
earthquake . 


' 


7 

The  first  investigations  in  the  U.S.S.R.  concerning  the  recording 
of  deep  reflected  waves  were  carried  out  by  G.  A.  Gamburtsev  and  others 
as  early  as  1939  (Beloussov  et  al.,  1962).  However  the  utilization  of 
such  waves  does  not  appear  to  have  been  realized  until  the  1950's  and 
the  development  of  deep  seismic  sounding  (DSS)  methods.  In  1956  detailed 
investigations  were  begun.  These  made  use  of  the  continuous  profiling 
method  with  a  large  number  of  shot  points  and  closely  spaced  seismometer 
groupings,  thus  enabling  simultaneous  study  of  the  sedimentary  strata  and 
the  deep  crust.  Such  deep  seismic  sounding  does  not  rely  on  any  one 
seismic  method  for  interpretation.  Wave  types  generally  used  are  the 
refracted  plus  sub-  and  super-critically  reflected  compressional  waves, 
but  some  efforts  are  made  to  employ  converted,  transverse  and  surface 
waves  as  well.  Attention  is  paid  to  both  amplitude  (dynamic)  and  velo¬ 
city  (kinematic)  properties  of  the  arrivals  in  order  to  achieve  the  best 
application  of  the  recorded  data.  A  survey  of  most  aspects  of  the  entire 
DSS  program  is  presented  in  a  book  recently  translated  from  the  Russian 
(Zverev,  1967).  In  this  monograph,  various  papers  indicate  Soviet  method¬ 
ology  up  to  about  1962.  Instrumentation  is  described,  data  is  presented, 
analysis  methods  are  shown  and  the  resulting  interpretations  are  given. 
This  collection  probably  provides  the  most  comprehensive  accumulation  of 
papers  concerning  the  Soviet  Union  crustal  study  program  presently  avail¬ 
able  in  English. 

In  an  earlier  paper,  Kosminskaya  and  Riznichenko  (1964)  provided 
a  summary  of  the  results  from  seismic  studies  in  various  parts  of  Russia. 
Their  review,  while  not  giving  the  details  that  were  possible  in  the 
book,  does  include  information  on  observational  methods,  on  the  use  of 
refracted  and  reflected  waves  and  on  the  crustal  cross  sections  obtained. 
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They  also  discuss  the  nature  of  crustal  stratification  in  a  section 
similar  to  a  paper  by  Kosminskaya  (1964).  A  very  brief  summary  of  the 
activities  in  the  DSS  field  by  all  East-European  countries  is  given  in 
a  report  prepared  by  Kosminskaya  (1965).  As  is  usual  in  the  Soviet  method 
of  crustal  investigations,  the  use  of  reflected  waves  is  dealt  with  as 
only  one  part  of  a  comprehensive  crustal  study  program. 

An  exception  to  this  generalization  is  an  article  by  Beloussov 
et  al.  (1962)  which  deals  specifically  with  the  recording  of  deep  reflec¬ 
ted  waves.  The  experiment  concerned  a  crustal  structure  study  using  a 
system  of  observations  that  ensured  a  continuous  correlation  of  sub- 
critically  reflected  energy.  It  is  worthwhile  to  mention  that  their 
method  required  a  grouping  of  sources  and  receivers  with  dimensions 
similar  to  those  the  present  research  found  necessary.  The  results  of 
their  investigations  indicated  a  number  of  reflecting  horizons  within 
the  lower  crust  and  even  in  the  upper  mantle.  But  their  detailed  sum¬ 
mary  of  the  reflections  shows  that  individual  phases  could  not  be  con¬ 
tinuously  correlated  over  more  than  about  ten  kilometers,  although  the 
'grouping'  with  depth  of  reflecting  horizons  was  concentrated  within  a 
few  kilometers.  Some  indication  of  fractures  in  the  deep  crust  was  given 
by  this  paper. 

In  the  first  article  of  the  monograph  previously  mentioned, 

Fursov  and  Yaroshevskaya  (1967)  include  a  section  on  the  recording  of 
sub-critically  reflected  waves  using  magnetic  tape  recorders.  They 
exhibit  a  number  of  records  with  reflected  energy  in  the  time  range  of 
12  to  16  seconds.  These  indicate  that  certain  amplifier  characteristics 
are  preferable  and  replaying  the  tapes  with  mixing  is  advantageous  for 
the  enhancement  of  reflected  signals.  Khalevin  et  al .  (1966)  also 
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point  out  that  mixing  improved  the  signals  from  deep  reflectors.  They 
show  no  examples  of  seismograms  but  give  a  crustal  cross  section  for  the 
Middle  Urals  based  on  DSS  data.  This  cross  section  depicts  the  crust  as 
being  horizontally  layered,  but  non-continuous ly ,  and  of  block  structure 
with  a  number  of  reflecting  horizons.  Zones  of  faulting  within  the  lower 
crust  and  on  the  Mohorovifcid  discontinuity  are  indicated.  Considerable 
structural  relief  is  suggested,  the  maximum  depth  variation  of  the  Moho 
lying  between  33  and  47  kilometers  over  the  extent  of  the  450  kilometer 
profile . 

While  geophysicists  in  the  Soviet  Union  have  been  developing  their 
DSS  program,  including  the  study  of  deep  reflected  waves,  German  geophysi¬ 
cists  have  been  the  most  active  in  using  such  waves  for  crustal  investi¬ 
gation.  In  cooperative  studies  among  geophysicists  from  oil  companies, 
governmental  agencies  and  university  groups,  thousands  of  reflection 
seismograms  with  coherent  energy  arriving  at  long  travel  times  have  been 
recorded  in  Germany.  One  of  the  first  reports  from  such  efforts  was 
that  of  Dohr  (1957)  in  which  more  than  250  deep  reflections  recorded  in 
the  course  of  routine  seismic  prospecting  were  observed.  Because  the 
deep  reflections  could  not  be  phase  correlated  over  more  than  a  few  kilo¬ 
meters,  Dohr  suggested  that  a  statistical  analysis  could  still  yield  a 
reasonable  interpretation.  Such  a  method  involved  the  construction  of 
histograms  in  which  the  number  of  observed  deep  reflections  within  a 
short  time  interval  (0.1  or  0.2  seconds)  from  seismograms  in  a  specific 
region  were  plotted  against  arrival  time.  Subsequent  to  this  work,  nearly 
all  German  interpretations  of  events  reflected  from  the  lower  crust  have 
followed  the  same  procedure.  A  question  might  be  raised  as  to  whether 
such  a  statistical  approach  was  absolutely  essential  in  these  cases  and 
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that  possibly  some  reflections  could  have  been  continuously  correlated 
over  considerable  horizontal  distances. 

An  exception  to  Dohr’s  approach  is  given  in  a  paper  by  Krey  et  al. 
(1961) .  In  an  investigation  using  seismic  reflection  measurements  to 
determine  whether  tectonics  of  the  sub-basement  could  be  correlated  with 
known  regions  of  ore  veins,  three  profiles  were  recorded  in  the  Siegerland 
district.  Reflections  from  the  immediate  sub-basement  were  generally 
weak  and  uncorrelatable ,  but  much  to  their  surprise,  reasonably  strong 
reflections  occurring  at  traveltimes  of  5  to  9  seconds  could  often  be 
well  correlated.  These  led  to  the  discovery  of  extensive  horizons  and 
of  large-sized  structural  features,  including  faults,  in  the  lower  crust 
(approximate  depths  of  10  to  25  kilometers).  Since  none  of  the  German 
review  papers  mention  the  methodology  of  recording,  other  than  it  being 
an  extension  of  routine  seismic  investigations,  it  is  worthwhile  to 
mention  the  experimental  techniques  revealed  by  Krey  et  al.  As  receiver 
arrays,  twelve  geophones  grouped  in  a  star- shaped  pattern  were  selected. 

The  source  generally  consisted  of  about  150  kilograms  of  explosive  per 
shot,  distributed  among  at  least  six  boreholes  of  about  15  meters  depth. 

The  experiment  was  unsuccessful  in  showing  a  relationship  between  the 
vein  features  and  tectonics  in  the  basement.  However  the  authors  state 
that  the  deep  structure  and  zones  of  fracturing  could  be  correlated  with 
the  regions  of  occurrence  of  the  ore  veins. 

In  a  general  review  of  crustal  structure  in  Western  Germany,  the 
German  Research  Group  for  Explosion  Seismology  (1964)  included  a  brief 
section  on  deep  reflection  investigations.  A  good  list  of  references  is 
given  and  one  15-kilometer  profile  is  presented.  Although  the  correlation 
of  one  phase  over  more  than  a  few  kilometers  was  not  possible,  the  profile 
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shows  the  accumulation  of  reflections  at  certain  levels.  Such  levels 
are  then  related  to  different  crustal  layers  interpreted  from  refrac¬ 
tion  seismology.  In  a  recent  paper  Dohr  and  Fuchs  (1967)  give  a  compre¬ 
hensive  review  of  deep  crustal  reflections  in  Germany.  By  displaying 
examples  of  seismograms,  they  indicate  that  it  is  preferable  and  advan¬ 
tageous  to  employ  groups  of  about  twelve  seismometers  with  higher  natural 
frequencies  than  refraction  instruments  in  order  to  successfully  record 
and  identify  deep  reflections.  Dohr  and  Fuchs  present  good  arguments  and 
evidence  for  the  existence  of  the  deep  reflections  in  light  of  criticism, 
which  is  enumerated,  against  the  reality  of  such  events.  Arguments  are 
also  made  in  favor  of  their  statistical  approach  to  the  interpretation 
and  some  questions  are  raised  as  to  the  physical  significance  of  their  re¬ 
sults.  In  the  course  of  interpretation,  it  is  assumed  that  multiple  reflec¬ 
tions  have  no  significance  at  the  long  traveltimes  of  their  observations. 

In  contrast  to  such  an  assumption,  our  investigations  have  shown  that  rever¬ 
berations  within  the  sedimentary  layers  can  produce  multiple  reflections 
which  have  significant  amplitude  within  one  second  of  the  onset  of  a  primary 
event  (Clowes  et  al.,  1968). 

The  question  of  the  physical  properties  of  the  deep  crustal 
reflectors  and  their  relation  to  the  observed  reflections  is  considered 
in  two  subsequent  papers  by  Fuchs  (1968,  1969).  In  the  more  recent 
paper,  he  concludes  from  the  German  data  that  levels  of  accumulated 
reflections  cannot  represent  first  order  discontinuities  nor  can  they 
represent  transition  zones  of  monotonically  increasing  velocity  with 
depth.  By  use  of  synthetic  seismogram  studies  and  other  arguments,  it 
is  shown  that  laminated  transition  zones  with  a  series  of  velocity 
reversals  and  increases  most  nearly  give  results  concordant  with  the 


observations . 
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In  order  to  supplement  and  confirm  the  findings  from  near- vertical 
incidence  profiling,  the  German  Research  Group  for  Explosion  Seismology 
(1966)  decided  to  record  a  wide-angle  profile  (sub-  and  super-critical 
reflections  from  0  to  150  kilometers)  about  a  common  subsurface  point. 
This  reference  describes  the  field  procedures  and  exhibits  some  record 
sections  with  reflections  from  the  Moho.  In  some  respects  this  work  is 
similar  to  the  initial  deep  reflection  program  at  the  University  of 
Alberta.  An  areal  pattern  of  shot  holes  was  used,  with  a  total  dynamite 
charge  ranging  from  45  to  200  kilograms.  Receiver  arrays  using  four 
vertically  recording  geophones,  with  group  spacings  of  about  100  meters, 
were  generally  laid  in  a  line.  In  order  to  distinguish  between  P  and 
SV  waves,  the  output  from  one  vertical-  and  one  horizontal-component 
geophone  at  the  same  location  was  recorded.  Meissner  (1966)  presented 
a  preliminary  interpretation  of  the  resultant  profile.  Four  different 
horizons,  corresponding  to  the  basement,  the  Conrad1  discontinuity,  a 
sub-Conrad  reflector  and  the  M  discontinuity  have  been  analysed  on  a 
T2  -  X2  plot.  The  strongest  amplitudes  were  for  reflections  from  the 
Moho.  From  a  consideration  of  irregular  traveltime  arrivals  which  did 
not  exactly  fit  the  expected  hyperbolic  curve,  Meissner  concludes  that 
velocity  gradients  over  about  two  kilometers,  as  opposed  to  sharp  dis¬ 
continuities,  explain  the  discrepancies  best.  In  a  later  paper, 

Meissner  (1967a)  considers  these  problems  in  a  more  complete  study.  He 
points  out  the  advantages  of  greater  amplitudes  for  wide-angle  measure¬ 
ments  than  for  near- vertical  ones,  especially  as  deeper  interfaces  are 
being  investigated.  This  is  especially  true  in  the  region  of  the 

1In  Europe,  the  Conrad  discontinuity  is  generally  regarded  as  the 

top  of  an  intermediate  layer  below  which  the  crust  is  more  basic . 
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critical  distance  where  the  writer  refers  to  the  large  amplitudes  expected 
for  "diving  waves"  (waves  for  which  the  turning  point  of  the  ray  is  within 
the  gradient  zone) .  Rather  strong  proof  is  given  for  the  existence  of 
these  diving  waves.  The  conclusions  from  this  study  indicate  that  the 
transition  zone  at  the  depth  of  the  M  discontinuity  is  comprised  of  a 
series  of  lamellas  with  material  of  stepwise  different  velocity  and  each 
less  than  about  150  meters  thick.  Meissner  points  out  that  temperatures 
at  this  depth  are  in  the  region  of  the  melting  points  of  most  components 
of  sialic  and  even  of  some  gabbroic  material.  Thus  he  postulates  that 
the  transition  from  sialic  to  basic  and  ultrabasic  matter  is  stepwise, 
interrupted  by  rhythmically  arranged  series  of  partial  melts  with  lower 
velocities.  For  a  more  comprehensive  study  of  the  interpretation  of  the 
wide-angle  measurements,  Meissner  (1967h,  in  German)  should  probably  be 
consulted. 

Since  the  review  by  Steinhart  and  Meyer  (1961),  further  instances 
of  the  recording  of  waves  reflected  from  within  the  crust  have  been 
reported  in  North  America.  Although  James  and  Steinhart  (1966)  have  up¬ 
dated  the  review,  some  results  are  worth  repeating  and  a  few  others 
should  be  mentioned.  Narans  et  al.  (1961)  carried  out  an  experiment  in 
Utah  to  determine  the  possibility  of  obtaining  normal  incidence  reflections 
from  the  lower  crust.  The  quality  of  the  seismograms  was  not  very  good 
so  a  statistical  approach,  similar  to  that  used  in  Germany,  was  applied. 
However  they  were  considering  only  ten  records  in  contrast  to  the  thousands 
examined  by  the  German  geophysicists.  From  the  histogram,  four  to  six 
possible  reflections  could  be  noted  and  two  of  these  were  at  times  which 
in  depth  would  correlate  well  with  known  refracting  horizons.  In  a  paper 
discussing  crustal  structure  in  New  Mexico  as  determined  from  the  GNOME 
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underground  nuclear  explosion,  Stewart  and  Pakiser  (1962)  suggest  that 
some  arrivals  can  be  attributed  to  sub-  and  super-critical  reflections 
from  within  the  crust  and  from  the  M  discontinuity.  This  identification 
is  based  on  reflection  traveltimes  calculated  from  the  crustal  model. 

By  extending  the  recording  time  in  routine  seismic  petroleum  prospecting, 
Robertson  (1963)  found  a  continuous  reflecting  horizon  in  the  sub-basement. 
He  argues  quite  effectively  that  the  arrivals  are  true  normal  incidence 
reflections  and  estimates  that  the  depth  of  the  horizon  varies  from  7.5 
to  14  kilometers  over  the  29-kilometer  length  of  the  profile.  Roller 
(1965)  reported  on  the  results  of  a  reversed  seismic  refraction  profile 
in  which  reflections  beyond  the  critical  angle  were  recorded  from  an 
intermediate  and  the  M  discontinuity.  In  an  interesting  study,  Sanford 
and  Long  (1965)  have  found  crustal  reflections  on  about  25  percent  of 
hundreds  of  high-magnification  seismograms  produced  from  microearthquakes 
near  Socorro,  New  Mexico.  The  best  interpretation  of  two  sharp  arrivals 
which  are  found  2.5  and  5.0  seconds  after  the  direct  S  wave  makes  them 
the  reflected  phases  SxP1  and  SxS1  from  a  discontinuity  in  the  crust  at 
a  depth  of  18  kilometers. 

Dix  (1965)  conceived  a  detailed  experiment  for  recording  near¬ 
vertical  incidence  reflections  in  a  region  of  the  Mojave  Desert  having 
a  very  thin  sedimentary  section.  Experimentally  the  procedure  used  three 
to  five  holes  per  shot  with  a  total  charge  of  100  to  300  kilograms. 
Fourteen-cycle  geophones  were  set  out  in  pairs  at  each  recording  station 
and  normal  reflection  prospecting  instrumentation  was  used.  To  ensure 

*The  nomenclature  for  these  phases  refers  to  transverse  waves  con¬ 
verted  to  longitudinal  waves  upon  reflection  at  a  depth  of  x  kilometers 
(SXP)  or  transverse  waves  reflected  at  a  depth  of  x  kilometers  (SXS). 
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the  arrivals  were  true  reflections,  some  cross-spread  recordings  were 
made.  Many  good  quality  reflections  were  observed  but  some  difficulties 
of  interpretation  are  mentioned.  One  strong  series  of  pulses  corresponds 
to  a  depth  of  25  kilometers,  the  depth  to  the  Moho  in  that  region. 

Perret  (1967)  reported  the  recording  of  vertical  incidence  deep  reflec¬ 
tions  on  four  special  high  sensitivity  vertical  accelerometers  in  a 
borehole  offset  just  166  meters  from  the  Salmon  underground  nuclear 
explosion.  The  data  were  recorded  on  magnetic  tape  and  digitized.  By 
time  shifting  and  compositing,  the  individual  reflections  are  shown 
quite  clearly.  In  addition  to  the  basement  contact,  interfaces  were 
determined  at  depths  of  13,  17.1,  S3  and  34.6  kilometers.  On  the  basis 
of  a  crustal  study  in  the  area,  the  last  depth  is  associated  with  the 
M  discontinuity  and  Perret  suggests  that  the  existence  of  the  two 
closely  spaced  deep  reflections  implies  a  transition  zone  at  the  crust- 
mantle  boundary. 

The  work  reviewed  in  the  last  two  paragraphs  has  provided  good 
evidence  for  the  existence  of  crustal  reflections  from  a  variety  of 
sources:  relatively  small  chemical  explosions,  nuclear  explosions  and 
near-earthquakes.  The  fact  that  reflections  are  recorded  from  such  dif¬ 
ferent  source  types  and  in  widely  separated  areas  strengthens  the  conten¬ 
tion  that  horizons  of  some  form  within  the  crust  are  able  to  reflect 
sufficient  coherent  energy. 

However  until  quite  recently,  the  University  of  Alberta  was  the 
only  North  American  group  systematically  using  the  seismic  reflection 
technique  for  the  investigation  of  crustal  structure  (Clowes  et  al„,  1968; 
Kanasewich  and  Clowes,  1968;  Kanasewich  et  al. ,  1969;  Clowes  and  Kanasewich, 
1969).  But  in  the  summer  of  1968,  Professor  R.A.  Phinney  of  Princeton 
University  began  a  program  which  would  use  reflected  waves  to  determine 
crustal  structure.  For  this  initial  work  one  crew  from  the  University 
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of  Alberta  joined  them  in  a  cooperative  experiment.  Reflections  from 
deep  within  the  crust  were  successfully  recorded  near  Riverton,  Wyoming. 
Realizing  the  success  of  the  technique  in  continental  areas,  the  group 
at  Princeton  undertook  a  similar  project  at  sea.  The  experiment  was 
successful  and  demonstrated  that  their  method  could  be  used  to  investi¬ 
gate  structure  at  least  to  depths  of  the  Moho.  Important  to  the  success 
of  the  project  was  the  utilization  of  a  receiver  array  with  low  noise 
level,  recording  of  data  on  FM  magnetic  tape  and  the  ability  for  subse¬ 
quent  digitization  which  allowed  the  application  of  signal  enhancement 
techniques,  particularly  beam  forming  (velocity)  filters.  Perkins 
(1969)  and  Phinney  and  Perkins  (1969)  have  reported  on  this  preliminary 
study.  It  is  worthwhile  mentioning  that  the  Seismology  Division  of  the 
United  Kingdom  Atomic  Energy  Authority  is  preparing  an  experimental 
program  to  record  reflection  arrivals  from  the  oceanic  crust  (Thirlaway, 
personal  communication).  As  well,  the  Bureau  of  Mineral  Resources  in 
Canberra,  Australia,  presently  has  a  group  actively  and  successfully 
working  on  deep  crustal  reflections  (Cleary,  personal  communication). 

The  detailed  seismic  studies  by  geophysicists  in  the  Soviet  Union, 
the  thousands  of  seismograms  examined  and  interpreted  by  researchers  in 
Germany,  the  success  of  the  program  at  the  University  of  Alberta  and  all 
the  other  preceding  documentations  suggest  that  reflections  from  within 
the  earth's  crust  can  be  obtained  and  used  to  determine  crustal  structure. 
Even  the  earlier  critics  of  such  work  (Steinhart  and  Meyer,  1961)  are 
accepting  this  contention,  but  with  some  reservations  (James  and  Steinhart, 
1966) .  However,  some  researchers  remain  unconvinced  that  reliable  sub- 
critical  reflections  are  recorded  from  the  M  discontinuity  (Knopoff ,  1969)  . 
Other  geophysicists  are  suggesting  that  the  method  of  reflection 
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seismology,  utilizing  both  normal  incidence  and  wide-angle  reflections, 
should  be  applied  to  a  much  greater  extent  in  the  future  in  order  to 
obtain  a  more  detailed  picture  of  the  anatomy  of  the  crust  (Francis, 

1969) . 

1 . 3  The  proj ect 

In  1964,  the  Geophysics  Division  of  the  Department  of  Physics 
started  an  experimental  program  to  investigate  the  possibility  of  obtain¬ 
ing  reflected  seismic  energy  from  boundaries  within  the  lower  crust  and 
from  the  Mohorovicid  discontinuity .  The  area  chosen  for  the  experiment 
was  an  east-west  line  in  southern  Alberta  where  the  general  crustal 
layering  was  known  from  the  results  of  reversed  refraction  profiles 
(Cumming  et  al.,  1962,  and  Maureau,  1964).  For  this  initial  work  a  modi¬ 
fication  of  the  expanding  reflection  profiling  technique  about  a  common 
depth  point  (Musgrave,  1962)  was  employed,  with  the  hope  of  obtaining 
average  vertical  velocities  to  any  reflecting  interfaces.  The  results 
of  this  preliminary  attempt  were  very  gratifying  and  were  reported  by 
Kanasewich  and  Cumming  (1965) . 

During  the  following  summer,  efforts  were  directed  toward 
extending  the  expanding  spread  profile  and  toward  obtaining  crustal 
reflections  in  regions  separated  by  as  much  as  120  kilometers.  While 
many  problems  were  encountered,  this  work  could  also  be  termed  successful. 
At  the  same  time  the  compilation  of  a  Bouguer  gravity  anomaly  map  for 
southern  Alberta  was  begun.  Subsequent  to  this  field  program,  some 
theoretical  seismogram  studies  of  transition  layers  were  initialized.  As 
well,  the  seismic  data  from  1964  had  been  digitized  thus  enabling  the 
application  of  signal  enhancement  techniques  by  digital  processing.  For 
this  preliminary  study,  compositing  of  traces  from  a  common  reflection 


o  ><  :  •  ,  Y-  o I  i 


via  aoiavrtqooD  91!?  ,  { <\Ql  rrl 


f  ■  i'  ' 


t?  oi  l  to  cOibutP.  mCTSOflISi  92  Ifil)  I  *9*109/1} 


18 

point  was  the  principal  technique  followed.  The  foregoing  paragraphs  are 
a  very  brief  summary  of  the  content  of  the  present  writer’s  M.Sc.  thesis 
(Clowes ,  1966) . 

Because  of  the  success  of  this  initial  experiment  and  because  some 
evidence,  from  both  seismic  and  gravity  anomalies,  had  been  obtained  for 
a  fault  in  the  lower  crust,  we  decided  to  make  an  intensive  study  of  this 
region  in  the  summer  of  1966.  The  method  by  which  this  was  accomplished 
was  continuous  seismic  reflection  profiling,  whereby  'continuous'  is 
meant  continuous  sub-surface  coverage  of  any  regions  which  might  reflect 
incident  energy.  The  resulting  profile  covered  approximately  40  kilo¬ 
meters,  beginning  at  the  NE  corner  of  Section  34,  Township  15,  Range  19 
W4M;  proceeding  north  but  angling  slightly  west  to  follow  roads  and  trails 
in  the  area;  finally  terminating  at  the  NE  comer  of  Section  31,  Township 
19,  Range  19  W4M.  To  this  line  was  appended  the  name  "Lomond  profile" 
and  it  is  so  indicated  in  Figure  1.1.  About  three  weeks  of  field  work 
were  necessary  to  obtain  the  data  along  the  Lomond  profile.  During  the 
same  summer,  further  gravity  station  locations  were  occupied  and  the 
data  from  these  were  added  to  the  previous  survey  and  the  general  compila¬ 
tion  of  gravity  values  kindly  supplied  by  the  Gravity  Division,  Dominion 
Observatory  of  Canada.  In  addition  a  start  was  made  on  covering  much  of 
southern  Alberta  with  a  network  of  ground  magnetometer  stations  to  supple¬ 
ment  the  gravity  and  seismic  data. 

The  reflections  obtained  from  the  Lomond  profile  were  of  good 
quality  and  showed  considerable  structural  relief  in  the  lower  crust  and 
probably  on  the  Mohorovi£id  discontinuity.  Thus  the  following  summer  it 
was  decided  to  extend  the  continuous  reflection  profile  northward  toward 
a  large  Bouguer  gravity  anomaly  which  suggested  further  possible  structure. 


.  ,  d  t  -  .•  M) 

f  Ins  zrnsnrxocpo  taui/ii  ;ir(?  So  ?  '  jdouz  orfj  io  o  u  d  a 

r  I  a  g  j  iiB  i  r?r  n.tTon  <  bt  ,do  •  M '  W 


joitioJ  of'it  gr.oLi  1  :  b  .  ;•  i  j*»o  o’  .  •>?.?  j;<5w 


'  fij  '  rf  .-1 


- 


19 


Figure  1.1  Location  map  for  the  seismic  profiles .  The 

heavy  dark  lines  are  the  two  continuous  reflection 
profiles .  The  dots  represent  locations  of  the 
recording  stations  for  the  refraction  survey . 
Bashed  lines  show  the  profiles  considered  for  the 
interpretation  of  the  refraction  data . 


>  P  "O*  -is  r.  '  •  •:/■>  m  ’ M 


LA  ?  j  • 


f  GREENBUSH 
■I  LAKE 
sf  50*  47'  N 
4  118*  21'  W  _Li 


L,iVf 
■  ■  111  s  * 


LAKE 


NEWELL 


““i  ^  | 


20 


To  avoid  a  difficult  crossing  of  the  Bow  River  it  was  necessary  to  offset 
this  second  line,  named  the  "Blackfoot  profile",  18  kilometers  to  the 
west.  Specifically  the  profile  started  at  the  NE  corner  of  Section  9, 
Township  19,  Range  21  W4M,  ran  north  for  14,5  kilometers  where  it  made  a 
3.2  kilometer  jog  to  the  west  and  then  proceeded  north  for  an  additional 
32  kilometers  terminating  at  the  NE  comer  of  Section  17,  Township  24, 

Range  21  W4M  (Figure  1.1).  Because  of  increased  efficiency  in  the  field 
procedure,  only  two  weeks  were  necessary  to  complete  this  set  of  data. 
Together  the  two  profiles  represent  about  90  kilometers  of  continuous 
sub-surface  reflection  coverage. 

Also  in  1967,  the  survey  with  ground  magnetometer  stations  was 
extended  in  southern  Alberta  and  into  southeastern  British  Columbia, 
thus  enabling  large-scale  magnetic  trends  to  show  up  more  clearly.  As 
well  it  was  found  possible  to  augment  the  gravity  compilation  with  data 
from  stations  which  had  been  occupied  in  southeastern  British  Columbia 
on  an  earlier  survey  by  the  University  of  Alberta  (Garland  et  al„,  1961). 

To  further  extend  the  methods  of  geophysical  analysis  of  the  area 
under  consideration,  a  cooperative  seismic  refraction  experiment  was 
made  with  personnel  from  the  Seismology  Division  of  the  Dominion  Observa¬ 
tory.  Four  4500-pound  charges  were  detonated  in  Greenbush  Lake,  British 
Columbia,  about  450  kilometers  west  of  the  reflection  profiles.  Some 
recording  stations  were  positioned  on  a  broadside  arc,  parallel  to,  but 
east  of  these  profiles,  at  locations  such  that  the  critically  refracted 
energy  would  be  incident  from  the  crustal  layering  directly  below  the 
profiles.  Other  stations  were  positioned  in  the  same  region  but  in-line 
so  that  apparent  velocity  values  for  the  refracted  waves  could  be  obtained. 
The  locations  of  the  recording  instruments  which  provided  useable  data 
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and  the  designations  of  the  profiles  are  also  shown  in  Figure  1.1. 

By  combining  the  results  of  these  various  methods  of  geophysical 
investigation,  it  proved  possible  to  confirm  interpretations  made  by 
any  one  method  and  at  the  same  time  extrapolate  detailed  structure 
obtained  from  the  reflection  data  over  a  more  widespread  area. 

1.4  The  technique 

Clowes  (1966)  described  in  detail  the  instrumentation  and  methods 
utilized  in  the  deep  crustal  reflection  program.  However,  certain  modi¬ 
fications  made  for  the  subsequent  field  projects  in  1966  and  1967  should 
be  mentioned  and  the  important  facets  of  the  method  can  withstand  repeti¬ 
tion. 

The  generally  good  quality  of  the  field  data  can  be  attributed 
principally  to  the  use  of  specially  designed  arrays  of  sources  and 
receivers.  Before  implementing  the  field  program,  theoretical  amplitude 
responses  of  various  combinations  of  hole  and  geophone  arrays  were  com¬ 
puted  using  a  method  adapted  from  antenna  theory  (Savit  et  al.,  1958). 
Figure  1.2  illustrates  the  linear  system  of  seismometer  and  hole  patterns 
which  were  generally  used,  together  with  their  normalized  amplitude 
responses.  The  combination  effectively  discriminates  against  incoming 
waves  with  apparent  wavelengths  small  relative  to  those  of  the  desired 
signals.  The  "relative  charge  size"  indicated  on  the  figure  was  deter¬ 
mined  by  a  binomial  distribution  which  approximates  the  desirable  spectral 
characteristics  of  a  Gaussian  function.  Typically  the  amount  of  explo¬ 
sives  planted  in  the  five  holes  was  2.3,  9.1,  13,  9.1  and  2.3  kilograms 
(5,  20,  30,  20  and  5  pounds),  respectively.  In  regions  where  the  energy 
coupling  of  the  explosion  with  the  ground  was  very  good,  these  quantities 
were  sometimes  reduced  by  one  half. 
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Figure  1.2  Normalized  amplitude  response  curves  for  a 
tapered  array  of  16  geophones >  6  holes  with 
differing  charge  sizes 3  and  the  geophone-hole 
combination .  The  arrow  at  indicates  the 
shortest  apparent  wavelength  of  the  desired 
signals  as  determined  from  seismograms . 


. 


23 


The  type  of  receiver  array  which  could  be  used  was  restricted  by 
the  number  of  geophones  per  take-out,  since  field  efficiency  had  to  be 
considered,  and  by  the  fact  that  only  pre-made  strings,  each  consisting 
of  four  geophones,  could  be  rented.  These  considerations  resulted  in  a 
'tapered'  array  of  16  geophones.  Twelve  seismometers  were  laid  along  a 
line  symmetric  about  the  main  cable  station  location  and  with  a  spacing 
of  12.2  meters  (40  feet)  between  detectors,  resulting  in  a  total  spread 
length  of  134  meters  (440  feet) .  An  additional  four  geophones  with  the 
same  spacing  were  laid  out  at  the  center  of  the  group  to  give  the  tapered 
array.  The  geophones  comprising  these  arrays  consisted  of  Electro-Tech 
EVS-4  7-1/2  Hz.  seismic  detectors.  At  this  point  it  should  be  reiterated 
that  Beloussov  et  al.  (1962)  successfully  used  a  system  of  arrays  much 
like  that  described  here.  They  regarded  the  grouping  of  seismic  sources 
and  receivers  as  a  necessity  for  the  isolation  of  sub-critically  reflected 
waves  and  the  suppression  of  interference  waves. 

For  the  project  of  continuous  profiling,  two  sets  of  truck-mounted 
instruments  each  recorded  the  ground  motion  from  six  arrays  of  geophones, 
the  12  resulting  traces  forming  one  half  of  the  normal  petroleum  explora¬ 
tion  'split-spread'.  Thus  two  shots  were  necessary  to  record  an  entire 
split-spread.  The  distance  between  arrays  was  set  at  292  meters  (960 
feet),  so  that  1.61  kilometers  (1  mile)  of  subsurface  coverage  was  recorded 
per  shot.  In  some  earlier  work,  twice  this  distance  was  used  but  it  made 
the  correlation  of  steeply  dipping  phases  across  one  record  too  difficult. 
The  data  from  the  arrays  were  simultaneously  recorded  on  photographic  film 
and  on  Precision  Instruments  6-channel  FM  tape  transports,  with  a  seventh 
channel  for  timing  purposes.  These  data  were  later  digitized  in  the 
laboratory  using  our  own  analog-to-digital  conversion  facilities.  In 
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addition  to  the  array  data,  three  Hall-Sears  HS-10  2  Hz.  geophones  were 
used  to  record  vertical  and  horizontal  components  of  earth  motion  on  the 
photographic  film.  In  general  filters  on  the  amplifiers  were  set  so  that 
energy  above  24  Hz.  was  attenuated  at  a  rate  of  24  db  per  octave.  The 
Texas  Instruments  recording  system  used  does  not  include  a  low-cut  filter 
so  the  low  frequency  attenuation  was  essentially  that  of  the  7-1/2  Hz. 
geophones.  No  automatic  gain  control  is  incorporated  into  the  instrumen¬ 
tation  so  that  data  between  0  and  4  to  6  seconds  generally  exceeded  the 
dynamic  range  of  the  amplifiers. 

The  experimental  method  just  described  proved  quite  successful 
for  recording  reflected  waves  from  the  lower  crust.  In  retrospect  it  was 
found  that  at  least  one  set  of  instruments  which  recorded  the  reflections 
over  a  wider  range  of  frequencies,  perhaps  by  a  factor  of  three  or  four, 
and  also  included  some  type  of  calibrated  automatic  gain  control  would 
have  been  extremely  valuable  for  the  subsequent  interpretation. 

1.5  The  field  observations 

More  than  100  seismograms  from  the  two  continuous  profiles  have 
been  recorded  photographically  and  on  FM  analog  tape.  The  quality  of 
these  field  observations  varied  from  very  good  to  very  bad,  but  generally 
they  were  of  acceptable  quality  (i.e. ,  reflected  energy  could  normally  be 
seen  on  the  photographic  monitors) .  Changing  shot-hole  conditions  with 
the  resultant  differences  in  energy  coupling  (each  set  of  boreholes  had 
to  be  used  at  least  twice)  and  variation  in  the  near-surface  unconsoli¬ 
dated  layers  on  which  the  geophones  were  planted  probably  contributed 
significantly  to  the  variation  in  quality  of  the  reflections.  One  further 
major  contribution  to  this  variation  was  likely  due  to  the  shape  of 
structures  being  investigated — concave  or  convex  reflectors  would  cause 
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focussing  or  spreading,  respectively,  of  any  reflected  energy.  Thus  the 
intensity  of  recorded  wavelets  could  vary  quite  rapidly  along  even  a  rela¬ 
tively  short  distance.  The  quality  of  some  of  the  seismograms  was  severely 
reduced  by  high  amplitude  60  Hz.  signals  whenever  it  was  necessary  to  lay 
a  string  of  geophones  near  a  commercial  power  line. 

In  Figure  1.3,  reproductions  of  some  of  the  originally  recorded 
seismograms  are  illustrated.  This  diagram  attempts  to  indicate  the 
variation  in  quality  of  the  recordings.  The  upper  three  sets  of  traces 
are  from  different  regions  on  the  Lomond  profile  and  show  examples  of 
good,  average  and  poor  quality  seismograms,  respectively.  The  lower 
record  is  an  example  of  a  good  quality  recording  from  the  Blackfoot  pro¬ 
file.  For  each  set  of  traces  illustrated,  the  upper  six  channels  are  the 
output  from  the  geophone  arrays  and  the  following  three  channels  represent 
ground  motion  obtained  from  single  2  Hz.  seismometers  recording  vertical, 
radial  and  transverse  motions,  respectively.  A  duplicate  of  these  nine 
traces  with  a  fourfold  reduction  of  gain  was  simultaneously  transcribed 
on  the  lower  half  of  each  record,  but  these  are  not  shown  in  the  figure 
except  on  the  poor  quality  seismogram:.  Since  the  instrumentation  does 
not  allow  automatic  gain  control,  such  channels  can  be  used  for  interpre¬ 
tation  if  the  attenuation  factors  were  set  too  low  and  excessive  amplitudes 
were  recorded  on  the  high  gain  traces.  Only  the  array  data  with  normal 
attenuation  were  recorded  on  magnetic  tape. 

As  evidenced  by  the  illustration,  it  was  usually  possible  to 
identify  some  reflections  on  the  photographic  seismograms  and  thus  make 
a  preliminary  interpretation.  Although  normal  moveout  is  insignificant 
at  the  traveltimes  shown,  in  many  cases  the  reflected  events  are  exhibit¬ 
ing  considerable  stepout,  in  both  directions  on  the  same  recording. 
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Figure  1.3 


Reproductions  of  deep  reflection  seismograms 


recorded  along  the  Lomond  and  Black  foot  profiles . 
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This  must  be  indicative  of  dipping  structures  in  the  deep  crust.  The  poor 
quality  data  set  shows  one  of  the  main  reasons  for  the  difficulty  in  ob¬ 
taining  reflection  information  suitable  for  i nterpretation — despite  the 
source-receiver  array,  large  amplitude  and  long  period  ground  motion  obscure 
the  reflected  energy  in  these  cases.  It  should  be  pointed  out  that  after 
the  seismogram  in  question  had  been  digitized  and  suitable  signal  enhance¬ 
ment  techniques  applied,  some  reflections  were  evident  on  the  resultant 
set  of  traces. 

As  previously  mentioned  the  special  refraction  profile  was  a  com¬ 
bined  effort  with  the  Seismology  Division  of  the  Dominion  Observatory. 

In  addition  to  providing  the  personnel  at  the  shot  point,  two  of  their 
crews  each  recorded  six  channels  of  seismic  data  on  either  Hall  Sears 
HS-10  or  Texas  Instruments  S-36  2  Hz.  seismometers,  these  being  separated 
by  0.5  kilometers.  Two  University  of  Alberta  crews  equipped  with  the 
same  instrumentation  as  for  the  reflection  studies,  except  for  the  use  of 
single  S-36  geophones  separated  by  0.29  kilometers,  each  recorded  twelve 
channels  of  seismic  information.  Three  of  these  were  the  vertical,  radial 
and  transverse  components  obtained  from  HS-10  seismic  detectors  located 
at  one  position.  Six  of  the  other  traces  were  also  recorded  on  the  FM 
analog  tape  transport.  One  additional  U.  of  A.  crew  registered  the 
vertical  and  horizontal  components  of  ground  velocity  as  produced  by  two 
Willmore  Mark  II  seismometers.  Twelve  of  the  seventeen  seismic  records 
were  considered  to  be  of  good  quality  for  interpretation.  Two  examples 
of  these,  one  from  a  Dominion  Observatory  (upper)  and  one  from  a  U.  of  A. 
crew  (lower)  are  illustrated  in  Figure  1.4.  On  this  diagram  the  initial 
arrival  of  the  head  waves  from  the  M  discontinuity  followed  by  a  later 
arrival  but  of  much  greater  amplitude  can  be  seen.  The  second  arrival 
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Figure  1.4  Reproductions  of  two  seismograms  recorded  in 
the  refraction  survey. 
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is  interpreted  to  be  a  head  wave  from  an  interface  within  the  Precambrian 
basement  complex.  This  small  refraction  project  provided  significant  cor¬ 
roborative  evidence  for  the  interpretation  made  from  the  reflection  seis¬ 
mograms  . 

The  following  chapters  will  show  that  valuable  and  useful  informa¬ 
tion  concerning  deep  crustal  structure  is  being  obtained  from  the  program 
of  geophysical  study  undertaken  in  southern  Alberta.  In  Chapter  2,  the 
application  of  digital  data  processing  techniques  to  the  reflection  seis¬ 
mograms  is  shown  to  have  considerable  importance  for  delineating  individual 
reflected  wave  trains  and  providing  an  improved  correlation  of  reflected 
events.  In  Chapter  3,  amplitude  spectra  of  the  observed  deep  reflections 
are  compared  with  the  spectra  from  synthetic  seismograms  in  order  to  give 
some  indication  of  elastic  properties  in  the  crustal  section.  Coincident 
with  this  study,  some  information  relating  to  the  physical  properties  of 
the  reflecting  horizons  (sharp  discontinuities  versus  transition  zones) 
is  provided.  Throughout  the  text,  the  feasibility  of  using  near- vertical 
incidence  reflected  energy  to  identify  horizons  in  the  lower  crust  is 
demonstrated.  In  contrast  to  the  assumptions  of  small  dips  and  planar 
interfaces  usually  incorporated  into  crustal  study  analyses.  Chapter  4 
provides  evidence  from  the  deep  reflection  seismograms  of  large-sized 
structural  relief  in  the  deep  crust  and  on  the  M  discontinuity,  structures 
which  rival  in  magnitude  the  largest  topographical  features  observed  at 


the  earth's  surface. 
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CHAPTER  2 

DIGITAL  PROCESSING 

2.1  Digitization  and  initial  processing 

All  the  seismograms  recorded  on  FM  analog  magnetic  tape  were  sub¬ 
sequently  digitized  by  the  technical  staff  of  the  geophysics  laboratory. 

A  block  diagram  of  the  system  is  illustrated  in  Figure  2.1.  The  digital 
tape  transports  indicated  on  the  illustration  were  Kennedy  incremental 
tape  recorders  having  seven  tracks  and  a  packing  density  of  200  charac¬ 
ters  per  inch.  In  order  to  provide  the  required  bandwidth  for  the  seis¬ 
mic  conversion,  the  original  analog  tape  drive  speed  of  30  inches  per 
second  was  reduced  prior  to  digitization  to  a  drive  speed  of  3-3/4  inches 
per  second,  a  factor  of  eight.  On  this  reduced  time  scale,  90  conversions 
per  second  were  made.  But  since  six  channels  were  being  sequentially 
digitized  (that  is,  N  points  for  each  channel  were  being  produced  as  in 
a  matrix  X—  where  i  varies  from  1  to  6,  as  j  goes  from  1  to  N) ,  this 
meant  that  each  individual  channel  was  being  sampled  at  a  rate  of  15 
conversions  per  second.  When  the  factor  of  eight  is  considered,  a  true 
conversion  rate  of  120  conversions  per  second  for  each  channel  is  realized. 
The  filters  used  prior  to  digitization  were  specially  designed  to  prevent 
aliasing  of  power  from  frequencies  greater  than  half  the  conversion  fre¬ 
quency,  that  is  the  Nyquist  frequency  for  the  digitized  data.  Figure  2.2 
shows  the  amplitude  response  of  the  filter  plotted  as  a  function  of 
frequency  as  recorded  in  the  field. 

The  digital  tapes  so  produced  in  the  laboratory  were  not  FORTRAN 
readable  on  the  IBM  360/67  computing  system  which  is  available  at  the 
university.  Consequently  a  machine  language  program  was  used  to  transfer 
the  data  onto  new  digital  tapes  and  in  a  format  such  that  the  information 
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Figure  2.1  Block  diagram  of  analog -to -digital  conversion 

system.  Blocks  of  data  are  read  alternately  in 
each  of  the  two  digital  transports  in  order  to 
provide  the  Inter-Record  Gap  required  for  the 
IBM  360/67  computing  system.  E.O.F.:  end  of 
frame ;  E.O.C.:  end  of  conversion. 
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F i gure  2 . 2 


Frequency  response  of  the  filters  used 
in  the  digitization  process. 
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would  be  FORTRAN  readable.  For  this  procedure,  the  two  tapes  containing 
the  digitized  data  from  six  channels  of  one  record  had  to  be  mounted 
simultaneously  on  two  tape  drives,  the  program  then  transferring  the  data 
onto  one  new  digital  tape.  All  necessary  checks  to  ensure  that  the 
seismic  information  was  being  transferred  as  originally  written  were 
included  in  the  program.  As  a  further  check,  parts  of  some  records  were 
plotted  and  compared  to  the  original  photographic  records.  In  this  manner 
all  the  digitized  seismograms  from  one  field  season  were  compiled  on  one 
FORTRAN  readable  tape.  A  duplicate  copy  of  each  tape  was  kept  as  a  safe¬ 
guard  against  possible  destruction  of  the  original. 

In  general,  a  total  of  either  thirty  or  fifty  seconds  of  two-way 
traveltime  starting  prior  to  the  shot  instant  and  continuing  from  that 
time  was  digitized.  The  timing  record  produced  in  the  digitization 
process  was  supposed  to  enable  the  determination  of  shot  instant  relative 
to  the  digitized  data.  However  more  often  than  not,  the  resultant  traces 
were  of  poor  quality  and  could  not  be  relied  upon.  It  became  necessary 
to  plot  some  of  the  data  from  each  record  which  was  on  digital  tape, 
compare  with  the  photographic  analog  records  on  which  timing  accuracy  was 
±2  milliseconds  and  in  this  manner  relate  which  point  corresponded  to 
zero  time.  This  procedure  proved  quite  reliable  and  it  is  estimated  that 
the  determination  of  shot  instant  is  generally  accurate  to  within  one 
digital  increment  or  ±8.3  milliseconds  and  certainly  not  more  than  twice 
this  amount. 

One  further  problem  was  illuminated  by  this  method  of  timing. 

The  process  of  digitization  took  place  over  periods  of  four  to  six 
weeks  in  two  different  years.  It  was  found  that  the  actual  sampling 
rate,  nominally  set  at  120  conversions  per  second,  varied  significantly 
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from  record  to  record  depending  on  when  the  actual  digitization  was  done. 
The  rate  varied  from  about  117  to  120  conversions  per  second.  Since 
some  of  the  analysis  techniques  which  were  to  be  applied  on  the  data  re¬ 
quired  a  constant  conversion  rate  and  for  the  sake  of  convenience  in 
data  handling ,  this  problem  was  corrected  as  discussed  in  a  later  para¬ 
graph.  It  was  a  tedious  process.  In  order  to  eliminate  such  problems 
in  the  future,  the  writer  suggests  that  the  shot  instant  be  recorded 
directly  on  one  channel  of  the  FM  tape  transport.  At  the  same  time, 
two  unique  and  identifiable  signals  (one  could  be  shot  instant) ,  care¬ 
fully  timed  so  that  they  are  exactly  separated  by  a  preset  value,  could 
be  generated  on  the  same  channel,  thus  allowing  for  a  quick  check  of  the 
digital  conversion  rate. 

Since  the  separation  of  the  source  and  receiver  arrays  usually 
varied  from  0  to  3.2  kilometers,  reflected  energy  did  not  arrive  at 
vertical  incidence,  even  for  plane  layers  at  the  base  of  the  crust.  For 
such  a  case  the  actual  variance  is  from  zero  to  less  than  three  degrees. 
The  difference  in  time  of  arrival  between  station  locations  due  to  non¬ 
vertical  incidence  is  called  normal  moveout  (NMO) .  Thus  the  initial 
processing  step  involved  the  removal  of  NMO,  effectively  reducing  the 
time  on  each  trace  to  that  for  vertically  travelling  waves.  While  doing 
this  it  was  considered  advantageous  to  select  a  specified  time  range  of 
interest,  usually  2  to  22  or  3  to  23  seconds  of  two-way  traveltime.  At 
the  same  time,  it  was  deemed  expedient  to  reduce  the  arithmetic  mean 
for  each  channel  to  zero: 
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Removal  of  the  NMO  was  actually  achieved  in  two  separate  steps. 
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First  a  computer  program  written  in  APL  (A  Programming  Language,  used  on 
remote  terminals  connected  to  the  central  processing  system)  was  utilized 
to  generate  the  time  differences  between  vertical  ray  paths  and  oblique 
incidence  ray  paths  for  all  times  in  the  range  of  interest  and  for  the 
entire  range  of  geophone-hole  offsets.  For  the  program  two  crustal 
models,  one  for  each  half  of  the  total  profile,  were  considered.  Velocity 
and  depth  information  for  the  sedimentary  layers  were  readily  available 
from  logs  of  wells  in  the  area.  Below  the  Paleozoic  sediments,  the  results 
from  reversed  refraction  profiles  and  the  preliminary  interpretation  of  the 
reflection  data  provided  the  necessary  information.  This  of  course  was 
a  case  of  requiring  a  preliminary  interpretation  before  subsequent  analyses 
could  provide  one  based  on  more  favorable  evidence.  However,  such  a  pro¬ 
cedure  simplified  the  computations  and  considering  the  low  sensitivity  of 
NMO  at  long  traveltimes  to  differences  in  crustal  structure,  even  for 
source-geophone  distances  of  a  few  kilometers,  the  method  was  deemed 
sufficiently  accurate.  Results  generated  by  this  APL  program  were  plotted 
as  a  form  of  normal  moveout  chart. 

Values  appropriate  for  each  channel  of  a  given  seismogram  were  read 
from  the  chart  and  fed  into  a  FORTRAN  language  computer  program  which 
performed  the  required  time  shifting  of  data  points.  Thus  a  new  tape 
containing  all  the  desired  records  from  one  field  season  and  consisting 
of  data  for  the  selected  time  range  with  d.c.  bias  and  NMO  removed  was 
generated.  In  another  suggestion  for  future  work,  it  is  felt  that  this 
two-step  procedure  could  be  combined  into  one  program,  thus  eliminating 
manual  plotting  of  the  NMO  chart,  reducing  the  number  of  data  input  cards 
required  and  facilitating  the  introduction  of  other  crustal  models. 
Alternatively,  if  sufficient  data  are  available,  corrections  for  NMO  can 
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be  applied  in  a  manner  analagous  to  such  determinations  in  the  seismic 
prospecting  field  (see  Brown,  1969). 

Previous ly,  mention  was  made  of  the  necessity  for  correcting  a 
variation  in  the  sampling  interval  from  one  seismic  record  to  another. 

Thus  the  data  with  NMO  removed  had  to  be  reprocessed.  Because  of  the 
small  differences  between  the  actual  digitizing  rate  and  the  specified 
standard  rate  of  120  conversions  per  second  and  to  make  the  calculations 
simple  and  efficient,  a  linear  interpolation  between  adjacent  data 
points  at  the  appropriate  time  was  employed  in  the  redigitization  pro¬ 
gram.  In  the  same  FORTRAN  program,  channels  showing  instrumental  dif¬ 
ficulties  were  zeroed,  those  with  reversed  polarity  were  corrected  and 
amplitudes  of  the  individual  channels  on  all  records  were  normalized 
within  a  specified  range.  For  this  normalization,  one  trace  considered 
to  have  a  good  average  amplitude  over  a  period  of  time  was  selected  as 
standard  and  all  other  traces  were  modified  to  give  a  similar  average 
amplitude  within  the  constraints  specified  in  the  program.  The  multi¬ 
plication  factor  for  each  trace  was  noted.  Such  a  process  is  not  ideal, 
but  with  the  variation  in  quality  and  amplitudes  of  the  original  record¬ 
ings,  it  proved  quite  effective  as  a  general  standardization  technique. 

Finally  a  complete  set  of  data,  redigitized  to  120  conversions 
per  second  with  NMO  removed  and  other  modifications  made,  was  obtained 
on  digital  tape.  For  the  20  seconds  of  recording  time  selected  for  each 
seismogram,  this  amounted  to  more  than  1.5  million  words  of  seismic 
information.  A  comparison  of  some  of  the  processed  seismograms,  repro¬ 
duced  in  analog  form  with  appended  time  scales  by  means  of  a  CALCOMP 
plotter,  was  made  with  the  original  photographic  analog  records.  For 
traces  near  the  shot,  where  NMO  would  be  a  minimum,  it  was  found  that 
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the  timing  of  identical  peaks  and  troughs  was  the  same  to  less  than  20 
milliseconds,  the  timing  accuracy  on  the  CALCOMP  plotted  data  being  about 
ten  milliseconds.  This  complete  initial  processing  procedure  was  tedious 
and  time-consuming,  but  absolutely  necessary  for  the  application  of 
signal  enhancement  techniques  and  the  efficient  handling  of  the  large 
amount  of  data. 

2.2  Autopower  spectra 

For  the  purpose  of  investigating  the  frequency  content  of  the 
reflection  seismograms,  a  FORTRAN  language  program  for  the  computation 
of  autopower  spectra  was  written.  The  actual  calculations  were  based  on 
a  fast  Fourier  transform  (F.F.T.)  algorithm  (Gentleman  and  Sande,  1966) 
for  which  the  basic  subroutine  programs  were  given  to  the  magnetotelluric 
analysis  group  of  the  Geophysics  Division  by  G.  Sande  and  provided  by 
them  to  the  writer. 

Basically  the  program  worked  in  the  following  manner.  After 
specifying  the  time  interval  of  one  seismic  trace  for  which  the  spectrum 
was  desired,  the  first  and  last  ten  percent  of  the  data  points  were 
tapered  with  a  cosine  bell  so  that  no  energy  from  sharp  steps,  which 
generate  power  at  all  frequencies,  would  be  introduced  into  the  spectrum. 
The  F.F.T.  algorithm  was  applied  to  produce  first  the  Fourier  coefficients. 
These  were  used  to  obtain  the  autocovariance  of  the  tapered  data  which 
was  then  modified  by  a  Parzen  lag  window.  The  principle  of  applying  a 
lag  window  to  the  autocovariance  is  to  provide  a  statistically  reliable 
power  spectrum  estimate,  that  is  one  for  which  the  variance  of  the  spec¬ 
tral  estimate  is  reduced.  One  characteristic  of  the  Parzen  window  is 
that  it  is  positive  definite  so  negative  power  estimates  should  never  be 
obtained.  In  addition,  only  a  fraction  of  all  the  autocovariance  values 
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are  required  for  the  calculations,  thus  increasing  the  efficiency  of  the 
computer  program  while  not  increasing  its  core  storage  requirements.  Af¬ 
ter  applying  the  window,  the  F.F.T.  subroutines  are  again  used  to  compute 
the  power  spectrum  estimate.  The  geophone  and  amplifier  responses  are 
removed  using  the  manufacturers  specifications  for  the  geophones  and  mea¬ 
sured  data  for  the  amplifiers.  Values  were  taken  at  the  appropriate  fre¬ 
quencies,  normalized,  squared  and  multiplied  with  the  power  spectrum  estimate 
to  yield  the  spectral  density  of  the  surface  particle  motion.  A  facility  is 
included  in  the  program  to  compute  the  autopower  spectrum  of  any  time  inter¬ 
val  on  individual  traces  or  provide  a  spectrum  averaged  over  several  channels 
of  one  record.  The  resultant  data  are  plotted  on  the  CALCOMP  plotter. 

Figure  2.3  illustrates  the  averaged  relative  autopower  spectrum 
of  four  different  seismograms,  these  being  chosen  as  representative  of 
the  general  spectral  characteristics  observed  from  an  analysis  of  all 
available  records.  The  two  spectra  on  each  plot  are  from  different 
traveltime  intervals,  as  indicated  to  the  right  of  the  curves,  and  cor¬ 
respond  to  increments  which  contain  good  reflections  from  the  R  discon¬ 
tinuity  (lower  curve),  an  intermediate  crustal  reflector,  and  from  the 
M  discontinuity  (upper  curve) .  In  order  to  facilitate  presentation,  the 
latter  has  been  shifted  by  adding  a  constant  to  each  estimate,  this  value 
also  being  noted  on  the  right.  The  ordinate  is  a  logarithmic  scale, 
but  since  the  autopower  amplitudes  vary  from  plot  to  plot  and  the  computer 
program  automatically  scales  these  to  fit  the  diagram  dimensions,  this 
scale  can  change  from  one  diagram  to  the  next. 

Record  07,  1966  (upper  right)  shows  a  feature  which  was  observed 
on  a  number  of  the  power  spectrum  plots— two  large  peaks,  one  usually 
centered  from  9  to  11  Hz.,  the  other  from  16  to  19  Hz.  This  same 
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Figure  2.3  Averaged  ccutopower  spectra  for  selected  time 
intervals  of  four  different  seismograms . 
Spectral  modifications  due  to  seismometer  and 
amplifier  characteristics  have  been  removed . 


RECORD  05  1966  RECORD  07  1966 

AVERAGED  SPECTRUM  FROM  6  TRRCES  AVERAGED  SPECTRUM  FROM  6  TRACES 
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tendency  is  shown  in  the  plot  for  Record  05,  1966  (upper  left)  but  the 
second  peak  is  much  subdued  relative  to  the  first.  As  well,  these  curves 
exhibit  the  lower  limit  of  frequency  range  (7  to  10  Hz.)  for  which 
spectral  peaks  were  observed.  In  contrast.  Record  25,  1967  (lower  right) 
shows  the  range  of  the  higher  limit  of  frequency  (12  to  17  Hz.)  observed 
for  the  principal  spectral  peak.  On  Record  55,  1967  (lower  left),  the 
peaks  do  not  appear  so  well  defined  but  their  central  frequency  is  about 
12  Hz.,  a  common  value  for  many  records.  As  a  consequence  of  examining 
these  and  similar  plots  for  all  available  seismograms,  it  was  found  that 
the  main  energy  of  the  reflected  arrivals  from  the  R  and  M  discontin¬ 
uities  lay  in  the  range  from  8  to  15  Hz.  Quite  a  number  of  recordings 
exhibited  a  doublet  peak  in  their  spectral  estimate.  It  was  also  noted 
that  the  central  frequencies  of  spectra  at  different  times  on  the  same 
seismogram  were  often  shifted  relative  to  one  another,  the  later 
reflection  usually  having  a  lower  central  frequency  (see  Record  05) . 
Another  aspect  of  the  general  character  of  the  plots  was  a  tendency 
for  the  higher  frequencies  to  have  a  lower  average  amplitude  than  the 
lower  frequencies.  These  observations  are  of  importance  to  the  study 
which  is  considered  in  the  next  chapter. 

2.3  The  zero-phase  band-pass  filter 

Since  the  autopower  spectra  showed  that  most  often  the  reflected 
energy  was  confined  to  a  relatively  narrow  frequency  band,  it  was  thought 
that  the  application  of  digital  band-pass  filters  would  be  a  suitable 
signal  enhancement  technique.  In  the  course  of  his  M.Sc.  thesis  research 
at  the  University  of  Alberta,  Alpaslan  (1968)  developed  a  FORTRAN  com¬ 
puter  program  for  a  zero-phase-shift  recursive  Butterworth  band-pass 
filter.  In  his  thesis,  a  good  summary  of  the  theory  is  given  and 
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references  are  provided.  The  Butterworth  filter  is  an  optimum  filter 
when  signal  and  noise  occur  in  separate  bands  of  the  frequency  domain 
and  does  not  exhibit  the  undesirable  Gibb's  effect,  a  high  frequency 
oscillation  superimposed  on  the  amplitude  spectrum  of  the  filter.  By 
utilizing  recursion  relations  for  the  computation  of  filter  coefficients 
and  the  application  of  these  to  the  data,  the  computer  program  is  made 
efficient  and  accurate.  A  desirable  characteristic  of  such  filters  when 
applied  to  real  data  is  that  they  be  zero  phase.  Otherwise,  a  differen¬ 
tial  lag  would  exist  between  different  frequency  components  thus  intro¬ 
ducing  distortion  into  the  filtered  data.  In  the  program  used,  the  zero- 
phase  aspect  was  introduced  in  a  simple  but  effective  manner.  The  output 
is  obtained  directly  from  a  product  of  the  z  transform  of  the  input  vector 
and  a  digital  filter  operator,  then  it  is  time  reversed.  The  time  re¬ 
versed  data  are  filtered  in  the  same  manner  as  the  original  input,  and 
finally  this  output  is  time  reversed  to  produce  the  desired  filtered 
trace  with  no  distortion.  Mathematically  it  can  be  shown  that  such  a 
process  does  indeed  realize  a  zero-phase  response,  the  only  significant 
difference  being  that  the  amplitude  response  of  the  filter  is  squared 
relative  to  that  of  a  single  pass  filter. 

An  eighth-order  band-pass  filter  of  the  type  just  described  was 
applied  to  a  number  of  digitized  seismograms.  Figure  2.4  shows  the 
original  recorded  data  (after  the  initial  processing  described  in  Section 
2.1)  for  four  records  and  the  same  seismograms  filtered  with  a  7  to  13  Hz. 
Butterworth  filter.  Considerable  improvement  is  rendered  by  the  process, 
but  an  important  disadvantage  was  noted.  Because  of  the  relatively 
narrow  frequency  specifications  of  this  filter  or  similar  ones,  the 
output  seismogram  appeared  quite  oscillatory  and  the  character  of  the 
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Figure  2.4  Blot  of  digitized  seismic  -reflection  records 

before  and  after  applying  a  7  to  13  Hz.  Butterworth 
digital  band-pass  filter.  The  illustration  is  on 
two  separate  pages. 

First  page:  seismograms  from  the  Lomond  profile . 
Second  page:  seismograms  from  the  Blackfoot  profile. 
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records  was  destroyed,  thus  making  phase  correlations  more  difficult. 

A  further  example  of  the  application  of  the  zero-phase  filter  to 
recorded  data  is  given  in  Figure  2.5.  The  upper  record  is  the  unfiltered 
seismogram.  Below  it,  the  result  of  applying  an  eighth-order  7  to  15  Hz. 
band-pass  filter  is  shown.  A  certain  amount  of  noise  reduction  is  evi¬ 
dent  especially  for  some  coherent  energy  about  14.4  seconds,  the  time 
corresponding  in  depth  to  the  M  discontinuity  as  derived  by  reversed  re¬ 
fraction  profiles.  The  autopower  amplitude  spectrum  of  Record  07,  1966, 
the  record  under  consideration,  showed  two  strong  peaks  centered  at 
about  10  and  17  Hz.  (see  Figure  2.3).  To  obtain  some  indication  as  to 
where  the  higher  frequency  energy  was  present  on  the  seismogram,  the 
original  record  was  processed  with  a  filter  having  cut-off  frequencies 
of  15  and  20  Hz.  The  third  set  of  traces  illustrates  the  results  which 
are  indicative  of  other  records  to  which  similar  filters  were  applied. 
Evidently  the  high  amplitude  sequence  at  about  11.6  seconds  contains 
most  of  the  energy  in  the  higher  frequency  range,  although  some  other 
envelopes  of  considerable  amplitude  can  be  observed.  Since  the  times 
about  which  the  envelopes  centre  are  coincident  with  the  observed  reflec¬ 
tion  times,  this  higher  frequency  must  be  intimately  related  to  the  re¬ 
flected  signal  and  should  not  be  isolated  or  attenuated.  As  will  be 
noted  later  in  Chapter  3,  considerable  information  relating  to  the  nature 
of  the  discontinuities  at  depth  is  present  in  the  spectral  density 
estimates.  Therefore  secondary  peaks  such  as  those  illustrated  in 
Figure  2.3  should  be  preserved  for  later  analysis,  not  filtered  out.  In 
the  light  of  our  experience,  it  can  be  stated  that  narrow  band-pass 
filtering  should  be  avoided  as  it  destroys  not  only  the  character  of  the 
reflected  impulse,  but  valuable  information  contained  therein. 
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Figure  2.5  An  example  of  different  filters  and  filtering 
techniques  applied  to  one  seismogram .  Upper: 
unfiltered  record;  second:  7  to  IS  Hz .  band-pass 
filter;  third:  IS  to  20  Hz.  band-pass  filter; 
fourth:  velocity  filter  (see  Section  2.4)  with 
a  pass  band  of  -35.2  to  +35.2  km/sec;  lower: 
the  velocity  filtered  output  passed  through  a 
7  to  IS  Hz.  band-pass  filter. 


fn_TER  butterworth  filter  butterworth  filter 


TIME  IN  SECONDS 


TIME  IN  SECONDS 

a. oo  a. 40  a. eo  9.20  9. so  10.00  10.40  io.ao  11.20  11. 60  12.00  12.40  i2.ao  13.20  13.60  14.00  14.40  14.80  is. 20  is. go  16.00 


j\/^yyyyyyyyy\/\/\/\A^\AAAAAAAAAAAAA^\AAAAA^---wwvNAAA/^AAAA^\^^~^^ — ~^aaa^wv\x 

- VWVWv — — ~W\ A/WWw< 


rvj  Q 

=1=  £ 


8.00  8.40  B.E 


TIME  (SECONDS) 

9.20  9.60  10.00  10.40  10. BO  11.20  11.60  12.00  12.40  12.80  13.20  13.C 


14.00  14.40  14.80  15.20  15.60 


LH 

cn 

o 


45 

Fuchs  (1969)  points  out  two  aspects  of  the  German  deep  reflection 
observations  which  should  be  mentioned  in  connection  with  these  results. 

He  refers  to  the  "band  character  of  the  reflected  signal"  and  the  shifting 
of  the  main  energy  from  one  phase  to  another  so  that  correlation  over 
more  than  a  few  kilometers  is  impossible.  Since  the  seismograms  in 
Germany  are  recorded  via  normal  seismic  reflection  prospecting  instrumenta¬ 
tion,  perhaps  the  filters  so  applied  cause  some  of  the  band  character. 

It  was  noted  earlier  that  even  the  application  of  a  low  frequency  band¬ 
pass  filter  from  7  to  15  Hz.  gave  an  oscillatory  appearance  to  the  resul¬ 
tant  seismogram.  From  the  results  of  processing  with  the  higher  frequency 
zero-phase  filter,  it  is  evident  that  correlation  of  one  phase  might 
prove  difficult  but  in  general  the  envelope  of  the  oscillations  can  be 
clearly  seen.  Thus  some  part  of  the  observed  results  in  Germany  may  be 
caused  by  the  recording  apparatus  rather  than  by  the  nature  of  the  re¬ 
flecting  horizons. 

2.4  The  velocity  filter 

In  the  process  of  analysing  the  photographically  recorded  seismo¬ 
grams,  it  was  noted  that  reflection  arrivals  with  varying  apparent  veloci¬ 
ties,  negative  and  positive,  could  be  seen  at  different  traveltimes  on 
one  record.  Quite  often  these  dipping  events  were  obscured  by  other  noise 
or  by  interference  with  one  another.  It  seemed  obvious  that  some  form  of 
filter  which  would  discriminate  against  all  reflected  energy  with  apparent 
velocities  outside  of  a  specified  range,  thus  enhancing  those  inside, 
could  prove  very  useful  in  distinguishing  among  various  dipping  events . 

Such  selective  filters  have  been  designed  and  applied  in  exploration 
seismic  prospecting  within  the  past  decade  and  have  been  called  by  various 
names:  velocity,  fan-pass,  pie-slice  or  beam- forming  filters  (for  example 
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see  Fail  and  Grau,  1963;  Embree  et  al,,  1963;  Wiggins,  1966;  Treitel  et 
al. ,  1967,  or  Sengbush  and  Foster,  1968), 

The  basic  idea  behind  such  filters  can  be  explained  quite  simply. 
By  designing  a  two-dimensional  time-space  operator  in  the  f-k  (frequency- 
wave  number  or  frequency-spatial  frequency)  domain,  it  is  desired  to  pass 
all  events  on  a  multichannel  seismogram  whose  apparent  velocities  fall 
within  the  range  -V  to  +V„  The  ideal  transfer  function  for  such  a  filter 
is 


Y  (  f ,  k  )  = 


+  1 


(2.1) 


I  0*  otherwise, 

where  k  is  the  spatial  frequency  related  to  the  wave  number  K  by  K  =  27rk, 
Figure  2,6  illustrates  the  ideal  filter  characteristic  of  equation  (2.1), 
Evaluation  of  the  two-dimensional  Fourier  transform  of  the  trans¬ 
fer  function  in  equation  (2,1)  will  produce  a  filter  operator  in  the  time- 
space  domain.  This  operator  represents  a  two-dimensional  array  of 
weighting  coefficients  which  must  be  convolved  with  the  data  points 
representing  traces  of  a  seismogram  in  order  to  generate  the  velocity 
filtered  output  data.  However,  Treitel  et  al ,  (1967)  have  derived  an 
algorithm  which  effectively  reduces  the  computational  labor  from  the 
original  formulation  of  the  problem.  They  then  make  use  of  a  recursive 
filter  (Shanks,  1967)  to  further  reduce  the  amount  of  calculations  re¬ 
quired.  Because  of  the  efficiency  of  this  method,  the  author  has  program¬ 
med  the  appropriate  equations  in  order  to  apply  this  form  of  filtering 
to  the  deep  reflection  data.  The  basic  procedure  can  be  briefly  described 
in  the  following  manner.  At  a  given  record  time,  the  digitized  data 
points  representing  a  specified  number  of  equally  spaced  input  traces 
are  appropriately  weighted,  time  shifted  and  summed.  This  weighted  sum 
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Figure  2„6 


Ideal  f-k  plane  transfer  function  of  a 
velocity  filter .  fp  and  kjj  are  the 
temporal  Nyquist  and  the  spatial  Nyquist 


frequencies respectively . 
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of  time-shifted  traces  is  then  convolved  with  the  fan  filter.  By  succes¬ 
sively  incrementing  along  the  vectors  of  data  points  and  repeating  the 
process,  one  velocity  filtered  output  trace  is  produced.  A  complete 
development  of  the  theory  is  given  by  Treitel  et  al.  (loc.  cit.).  The 
Appendix  contains  a  listing  of  the  computer  program  to  implement  such  a 
velocity  filtering  process,  with  adaptability  for  application  of  three 
different  filter  characteristics  to  each  set  of  input  data  and  for  CALCOMP 
plotting  of  all  results. 

Previously  it  has  been  mentioned  that  six  seismic  channels  from 
each  of  two  truck-mounted  instruments  were  recorded  for  each  shot,  these 
twelve  traces  then  yielding  one  mile  of  sub-surface  coverage.  While  the 
program  allowed  a  velocity  filter  with  a  range  (even  number)  of  traces, 
it  was  decided  that  the  application  of  an  eight-trace  fan-pass  filter 
would  be  adequate.  This  generated  five  filtered  output  traces  for  every 
twelve  input  channels. 

An  example  of  the  effectiveness  of  this  two-dimensional  operator 
in  passing  events  with  apparent  velocities  in  the  specified  range  is  given 
in  Figure  2.7.  Seven  different  dipping  events  with  apparent  velocities 
ranging  from  +11.7  km/sec  through  an  infinite  apparent  velocity  and  down 
to  -11.7  km/sec  (stepouts  of  +3  to  -3  sampling  intervals  per  trace)  are 
synthesized  in  the  input  channels  which  form  the  upper  part  of  the  figure. 
The  three  lower  sets  of  five  traces  each  are  the  result  of  processing  the 
input  data  with  the  three  different  pass  bands  indicated  on  the  diagram. 

It  is  evident  from  this  illustration  that  similar  dipping  events  with 
apparent  velocities  inside  the  specified  range  are  passed  with  no  modifi¬ 
cation  and  events  whose  apparent  velocities  fall  on  the  cutoff  velocities 
have  their  amplitudes  attenuated  by  a  factor  of  one  half.  All  reflec¬ 
tions  whose  apparent  velocities  lie  outside  the  range  of  the  fan-pass 
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Figure  2.1  Application  of  the  velocity  filter  to  synthetic 
data .  Upper:  twelve  synthesized  traces  with 
coherent  events  having  a  range  of  apparent  velo¬ 
cities  from  +11.7  through  0  to  -11.7  km/sec. 

Lower:  the  results  of  three  different  pass  hands 
(+35.2  to  +11.7  km/seCf  -35.2  to  +35.2  km/seCj  and 
-11.7  to  -35.2  km/sec3  respectively)  applied  to 
the  synthetic  data .  An  eight-trace  filter  was 
used o 
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filter  are  severely  attenuated.  It  should  be  mentioned  that  the  resolu¬ 
tion  of  the  velocity  filtering  method  is  limited  by  the  sampling  incre¬ 
ment  in  the  time  domain  and  by  the  distance  between  adjacent  station 
locations  in  the  space  domain.  For  the  present  research,  the  resolution 
limits  were  found  to  be  acceptable,  but  if  more  precise  determination  of 
dips  of  specific  events  were  desired,  one  or  both  of  the  above  parameters 
would  have  to  be  reduced. 

An  example  of  the  application  of  the  same  filters  (+11.7  to  -11.7 
km/sec)  to  a  set  of  real  data  is  given  in  Figure  2.8.  In  the  upper  part 
of  the  illustration,  the  original  twelve  traces  (after  initial  processing) 
which  have  a  spatial  separation  of  0.29  kilometers  are  plotted.  The 
velocity  filtered  seismograms  are  drawn  below,  all  amplitudes  remaining 
constant  for  each  plot.  The  good  quality  reflections  labelled  R  and  M 
were  extracted  from  the  original  data  by  a  filter  whose  apparent  velocity 
pass  band  was  from  -11.7  km/sec  to  -35.2  km/sec.  Two  other  possible  re¬ 
flections  with  less  amplitude  can  be  found  at  9.2  and  12.4  seconds. 

These  are  not  nearly  so  clearly  delineated  by  the  other  filters.  However 
a  reflected  event  with  a  different  dip  at  about  8.6  seconds  is  shown  on 
the  lower  two  sets  of  output  traces „  Other  coherent  phases  of  lesser 
amplitude  are  seen  throughout  the  filtered  sections.  Yet  one  must  be 
cautious  in  interpreting  such  output  since  random  noise  could  produce  a 
coherent  line-up  of  energy  which  would  be  accentuated  by  the  filter.  In 
general,  continuity  over  a  few  adjacent  sets  of  filtered  data  was  neces¬ 
sary  before  such  events  would  be  interpreted  as  possible  reflections. 

An  example  of  the  fan-pass  filter  applied  to  some  good  quality 
data  was  included  in  the  fourth  set  of  traces  displayed  in  Figure  2.5. 

For  this  case  the  pass  band  of  apparent  velocities  was  from  -11.7  to 
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Figure  2.8  Application  of  the  velocity  filter  to  real  data „ 
First:  the  unfiltered  input  traces;  second:  pass 
band  from  -11,7  to  -35,2  km/sec;  third:  pass  band 
from  -36 o 2  to  +35,2  km/sec;  fourth:  pass  band 
from  +35 o 2  to  +11,7  km/sec „ 


IS  35.?  K/s  VI  IS -3S.2  «/s  VI  IS  -11.7  n/s  VELOCITY  RECORDS  Ml R. 42R  1966 

FILTERED 

IS  11.7  K/S  VJ^  15  35. 2  K/S  V2„  IS-3S.2K/5  PCTP  DISTANCE  BETWEEN  TRACES  IS  .29  NM. 


8 


10 


R 


12 


M? 


14 


16 


VwvAotyWMA/"-7^'* 
_  vWvV(M^WWwv 
■wWvvV MAM/vWVw^ 


52 


-35.2  kilometers  per  second.  It  was  initially  thought  that  band-pass 
filtering  of  such  output  would  provide  a  very  good  seismogram.  To  this 
end  a  7  to  15  Hz.  Butterworth  filter  was  applied  to  the  velocity  filtered 
data.  As  the  last  set  of  traces  on  Figure  2.5  shows,  the  problem  con¬ 
cerning  loss  of  reflection  character  is  still  significant.  It  was  con¬ 
cluded  that  little  improvement  of  reflection  identification,  and  possibly 
a  loss  of  information,  would  be  rendered  by  applying  band-pass  filters 
to  the  fan-pass  filtered  output. 

Figure  2.9  is  an  illustration  of  the  tremendous  improvements 
which  can  be  attained  by  application  of  processing  techniques  described 
in  this  chapter.  The  upper  records  are  reproductions  of  two  photograph¬ 
ically  recorded  seismograms  which  together  form  one  record  as  required 
for  the  velocity  filtering  process.  They  are  poor  quality  seismograms 
from  which  no  reliable  information  could  be  derived,  even  when  compared 
with  adjacent  seismograms.  In  the  process  of  digitization  of  the  cor¬ 
responding  analog  tapes,  the  60  Hz.  noise  which  overrides  some  of  the 
traces  is  attenuated.  The  initial  processing  procedures  described 
earlier  include  a  gain  factor  applied  to  the  amplitude  of  the  traces. 
After  processing  by  a  fan-pass  filter  in  the  apparent  velocity  range  of 
35.2  to  11.7  km/sec  and  the  application  of  further  amplitude  gain,  the 
five  traces  of  velocity  filtered  output  for  the  same  data  as  the  upper 
seismograms  is  illustrated  in  the  lower  part  of  the  figure.  The  improve¬ 
ment  is  self-evident.  Reflection  events  with  significant  amplitude  are 
seen  throughout  the  plot  and  most  of  these  could  be  correlated  with 
similar  phases  on  adjacent  traces.  Needless  to  say,  the  improvement  of 
all  records  by  the  application  of  signal  enhancement  techniques,  and 
particularly  velocity  filtering,  was  not  so  spectacular  as  that  of  the 


. 
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Figure  2.9  An  example  of  signal  enhancement  by  digital 
data  processing  techniques .  Upper:  repro¬ 
ductions  of  two  field  seismograms  obtained 
from  one  shot .  Lower:  the  resultant  seismo¬ 
gram  after  digitization,  amplitude  gain  and 
application  of  a  velocity  filter  with  a  pass 
band  of  +35.2  to  +11.7  km/sec. 
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present  example.  However,  such  an  illustration  emphasizes  the  consider¬ 
able  importance  of  data  processing  methods  for  improvement  of  the  quality 
of  deep  reflection  seismograms. 
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CHAPTER  3 

CRUSTAL  ATTENUATION  AND  THE  NATURE  OF 
THE  REFLECTING  HORIZONS 

3,1  The  synthetic  seismogram  with  spatial  attenuation 

An  important  parameter  of  the  solid  earth,  and  one  to  which  much 
attention  has  recently  been  given,  is  the  attenuation  of  seismic  energy. 
If  this  fundamental  property  were  absent,  the  energy  of  past  earth¬ 
quakes  would  still  be  reverberating  throughout  the  earth  today.  Never¬ 
theless,  in  some  studies  it  is  found  that  agreement  between  theoretical 
results  and  observational  data  is  sufficient  to  a  first-order  approxima¬ 
tion  without  including  an  attenuation  factor  in  the  theory.  With  the 
recent  availability  of  broad-band  systems  capable  of  recording  a  large 
dynamic  range,  it  is  possible  to  obtain  data  by  which  such  a  factor  may 
be  studied.  As  a  measure  of  the  anelasticity  of  a  medium,  Knopoff  (1956) 
first  introduced  Q,  the  specific  attenuation  factor  or  dimensionless 
quality  factor,  into  the  field  of  seismology.  The  quantity,  Q,  repre¬ 
sents  a  reduction  to  a  dimensionless  form  of  the  more  usual  measures  of 
attenuation.  All  of  these  definitions  are  related  to  an  expression  in 
electrical  circuit  theory  which  gives  the  relative  efficiency  of  energy 
storage  at  resonance  (Ryder,  1959) : 


where  AE  is  the  energy  dissipated  per  cycle  and  E  is  the  maximum  energy 
stored  per  cycle.  To  relate  more  closely  to  seismology,  AE  can  be 
considered  as  the  amount  of  energy  dissipated  per  cycle  of  a  harmonic 
excitation  in  a  certain  volume  and  E  is  the  maximum  elastic  energy  in 
the  system  in  the  same  volume  (Knopoff,  1964) .  The  various  means  of 
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defining  Q  are  related  to  the  type  of  observation  being  used — for  example, 
standing  waves  versus  progressive  waves.  Knopoff  (loc.  cit.)  discusses 
the  related  definitions  which  are  identical  for  homogeneous  and  isotropic 
media. 

In  his  earlier  article,  Knopoff  (1956)  included  Q  in  the  differ¬ 
ential  equation  of  motion  of  a  medium  possessing  solid  friction.  Simpli¬ 
fied  to  the  case  of  plane  compressional  waves  propagating  in  the 
z-direction,  the  equation  of  motion  becomes 

[  i  +-L-  i-  ]  =  1-  if**  ,  (3.i) 

|  a)  |  Q  9t  9z2  c2  9t 2 

where  go  is  the  angular  frequency  whose  absolute  value  is  taken  to  guaran¬ 
tee  absorption  of  both  positive  and  negative  frequency  components  of  the 
Fourier  spectrum,  w  is  the  displacement  of  a  particle  and  c  is  the  wave 
velocity  had  the  losses  been  assumed  to  be  absent, 


c  =  (  L±J»  )  1/2  . 

Here  A  and  y  are  the  Lam6  elastic  constants  and  p  is  the  density.  If 
it  is  assumed  that  w  =  We^*3,  equation  (3.1)  becomes 

r  ,  iw  -i  d2W  go2  IlT  A 

L  1  +  -  J  -  +  —  W  =  0  , 

|go|Q  dz2  c2 

which  has  the  solution 

-1/2 

w  =  exp  [  iwt  -  i  —  (  1  +  ■  )  z  ]  . 

c  Wq 

If  Q  >>  1,  the  factor  in  parentheses  can  be  expanded  by  the  binomial 
expansion  (omitting  go/  |  go  | )  : 
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w 


exp  {  iwt 


io)Z  (jOZ 

c  "  2cQ 


.  i<25.  [  .  -i_  +  o  (  I  )  ] 

c  8Q2  Q4  (3.2) 

-22-  0  (  —  )  }  . 

2cQ  Q2 

Equation  (3.2)  shows  that  the  attenuation  is  approximately  expressed  by 
e^cQ  and  if  the  spatial  attenuation  factor  is  denoted  in  the  usual 

"  OtZ 

manner  by  e  ,  then  the  value  of  a  is  given  by 

°  =  2^  •  (3'3) 

In  this  manner,  the  specific  attenuation  factor,  Q,  is  related  to  the 
spatial  attenuation  coefficient,  a,  for  travelling  waves.  Q  can  also  be 
related  to  other  parameters  for  different  observational  data.  In  fact, 
the  determination  of  the  Q- structure  of  the  earth  is  an  important  study 
to  provide  an  additional  quantity  from  which  the  physical  and  chemical 
state  of  the  earth's  interior  may  be  inferred.  Sato  (1967)  has  given  a 
comprehensive  review  of  Q.  and  its  determination  from  observational  data. 

In  contrast  to  the  large  amount  of  attention  which  has  recently 
been  given  to  fundamental  studies  of  Q  in  the  earth's  interior,  there  has 
been  relatively  little  application  of  the  analysis  in  generating  synthetic 
reflection  seismograms.  The  synthesis  of  seismograms  is  a  technique  which 
has  only  been  available  since  the  1950 's  when  continuous  velocity  logs 
became  available  (see  for  example,  Vogel,  1952  or  Summers  and  Broding, 
1952).  Subsequently,  Peterson  et  al.  (1955)  described  an  analog  computer 
system  for  the  generation  of  synthetic  seismograms.  They  were  followed 
by  other  researchers  who  described  various  methods  by  which  seismogram 
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synthesis  could  be  achieved,  especially  by  use  of  digital  computers 
(Berryman  et  al.,  1958;  Diirschner,  1958;  Wuenschel,  1960  and  Bois  et  al . , 
1960).  Common  to  all  these  methods  was  an  assumption  of  perfect 
elasticity.  Baranov  and  Kunetz  (1960)  proceeded  one  step  further  and 
provided  means  by  which  depth- independent  attenuation  could  be  included. 
Trorey  (1962)  points  out  the  limitations  such  assumptions  impose  on  the 
theoretical  seismograms  and  proposed  a  method  by  which  frequency  and 
depth-dependent  absorption  could  be  programmed  into  the  computations. 

His  time-domain  solution  is  a  modification  of  the  ray-tracing  scheme  of 
Baranov  and  Kunetz  (loc,  cit.)  in  which  the  seismic  section  is  divided 
into  many  layers  of  equal  time  thickness. 

For  the  study  of  synthetic  seismograms  of  the  deep  reflection 
data,  it  is  important  to  be  able  to  investigate  the  effects  of  transition 
zones,  which  may  be  described  as  layers  of  the  earth  in  which  the  velocity 
increases  or  decreases  linearly  with  depth.  Berryman  et  al.  (1958) 
developed  an  iterative  formula  for  calculation  of  the  reflection 
coefficient  at  the  surface  for  a  model  section  comprised  of  horizontal 
transition  layers.  Plane  waves  at  normal  incidence  to  the  layers  were 
assumed  and  the  density  was  considered  to  be  a  constant.  In  the  course 
of  research  for  the  present  author's  M.Sc.  thesis,  the  equations  deduced 
by  Berryman  et  al .  were  programmed  to  calculate  the  transfer  function 
of  the  layered  system.  Since  the  crustal  section  is  known  to  attenuate 
seismic  energy  and  from  some  observations  made  when  the  autopower  spectra 
of  the  digitized  seismograms  were  investigated  (Section  2.2.),  it  was 
concluded  that  synthetic  seismogram  studies  incorporating  a  depth- 
dependent  attenuation  factor  were  necessary.  Consequently,  the  solution 
given  by  Berryman  et  al .  (1958)  has  been  rederived  in  order  that 
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attenuation  can  be  included  in  the  computations. 

The  equation  of  motion  for  the  displacement,  u  =  (ui ,  U2,  U3;t), 
of  a  point  in  an  unbounded  elastic  solid  due  to  the  transmission  of  a 
time-dependent  stress,  p^ . ,  at  any  instant,  t,  can  be  written  as 


92u. 


at' 


9p .  . 

13 

9x . 

3 


(3.4) 


where  p  is  a  constant  density  and  the  summation  convention  is  used. 


The  stress  is  obtained  from 


Pij  =  A  0  6 .  .  +  2  y  e .  . 

13  13 


where  0  is  the  cubical  dilatation. 


0  = 


9u. 

1 

977  * 

1 


i)  . 


and  e. .  is  the  strain  tensor, 

13 

.  9u.  9u. 

e  =  I  (  __J.  + 
ii  2  1  3x.  3x. 

1  3 

In  addition,  it  is  assumed  that  A  =  A  (xj,  x2,  X3)  and  y  =  y  (xi,  x2,  X3) 
are  spatially  dependent.  Substituting  into  (3.4)  and  neglecting  the  second 
and  higher  derivatives  of  the  Lam6  functions  gives  the  equation 
pU  =  (A+y) V0  +  0VA  +  yV2u  +  (Vy.V)u  +  V^(Vy.u)  ,  (3.5) 

When  studying  reflection  data  at  near-vertical  incidence,  a  use¬ 
ful  approximation  is  an  assumption  that  plane  waves  propagate  only  in 
the  z-direction.  In  this  case,  the  only  displacement  is  U3  =  w  and  only 
derivatives  with  respect  to  X3  =  z  are  non- zero  and  need  to  be  considered. 
By  again  neglecting  the  second  derivative  of  y,  equation  (3.5)  becomes 
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Let  w(z,t)  =  W(z)ei^t 

Substituting  into  equation  (3.6)  gives 

-M2  W(z)  =  4-  [  M*)  ]  dw 

dz  p  dz 


(3.7) 


+  r  X(z)  4-  2p(z)  d2W 
L  p  J  dz7 

where  the  total  differential  is  written  because  W,  X  and  y  are  functions 
of  z  only.  In  order  to  introduce  attenuation  into  the  problem,  X  and  y 
are  formally  expressed  as  complex  quantities  (Sato,  1967): 


X  =  XR  +  iX1 
y  =  y  +  l  yx 

As  the  density  is  constant,  all  the  attenuation  is  associated  with  these 
complex  elastic  modulii  and  the  square  of  a  complex  velocity,  v,  can  thus 
be  written: 


v 


2  _ 


XR  +  2yR  +  i(xl+2yl) 


(3.8) 


Substituting  this  into  equation  (3.7)  gives 


_a)2  W  ”  dz  ^  ^ 


dW 

dz 


+  v2(z) 


d2W 

dz2 


_d_ 

dz 


[  v2(z) 


dW  -| 
dz  ■* 


Upon  rewriting,  the  equation  becomes 


■$£  (  v2  §  )  +  “2w  "  0  '  (3'9) 
Consider  a  layered  medium  comprised  of  N+l  transition  layers.  The 
complex  velocity  in  the  j-1-*1  layer  may  be  written  as  a  linear  function  of 
depth : 


v(z) 


vj +  bj  <z-zP 

VjR  +  i  Vjl  +  bj 


Zj  ^  z  ^  zj+l 

(z-Zj). 


(3.10) 
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For  the  zeroth  layer,  there  is  no  upper  boundary  so  Zq  is  arbitrarily 
set  equal  to  zero.  In  this  model  the  initial  velocity  is  a  constant, 
Vq  =  Vi .  The  layer  is  an  infinite  half-space  of  velocity  vN.  Thus 
for  a  given  layer,  the  real  part  of  the  velocity,  representing  the 
travelling  wave,  has  a  linear  gradient  and  the  imaginary  part  of  the 
velocity,  representing  the  attenuation,  is  constant. 

Writing  equation  (3.9)  for  the  jth  layer  yields 


([(v-^+ivj1)  +  bj  (z-z.)]2^)  J  +  “2wj  •  °- 


(3.11) 


Equation  (3.11)  is  similar  to  equation  (5)  of  Berryman  et  al .  (1958) 


except  that  Vj  is  complex: 

v 


3  =  vjR  +  ivjT 


(3.12) 


The  development  in  their  paper  may  then  be  followed  step-by-step  to 


produce  the  reflection  coefficient  Rj  ^  for  the  (j-l)th  layer  in  terms 


of  Rj  . 


R 


l-R-j 

^  (l+R-  ) 


.1 


j-l  =  b.j-lgj-l  +  (bP  -  bJ6?  1+RJ  oxpC-B.  !  ln[Vj  -+..  -1iV3  ] 

1-R.  * 

bj-l6j_l  -  (bj  -bj.x)  +  bj  ( 1+Rj) 


vj-l  +  1V 


.  I 
"j-l 


(3.13) 


where 


6- 


[  1  -  4w2/bj 2  ]1/2  . 


For  the  layer  there  is  only  a  transmitted  wave,  so  Rjj  =  0.  On  a 

digital  computer,  equation  (3.13)  can  be  iterated  from  R^j  to  produce  RQ 
which  represents  the  ratio  of  the  amplitude  of  the  reflected  wave  at  the 
’surface’  to  that  of  the  incident  wave.  When  determined  for  all  frequen¬ 
cies,  this  constitutes  the  transfer  function  for  reflected  waves  at 
vertical  incidence  for  the  entire  layered  system.  All  multiple  reflec¬ 
tions  are  inherent  in  the  computed  result. 
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In  the  practical  application  of  equation  (3.13),  some  questions 
might  be  raised  as  to  what  constitutes  the  complex  velocity.  The 
available  data  are  continuous  velocity  logs  for  the  sedimentary  section, 
supplemented  by  velocity  and  depth  determinations  in  the  deeper  crust. 
These  values  give  the  phase  velocities  for  propagation  of  seismic  waves . 
Equation  (3.13)  requires  the  complex  velocities  so  a  relation  between 
the  specific  attenuation  factor,  (^,  and  the  real  and  imaginary  parts 
of  the  velocity  is  necessary.  In  general,  the  wave  equation  has  a 
solution  of  the  form  exp  [ico  (z/v  -t)]  where  v  =  vr  +  i  Vj.  (For  the 
sake  of  the  following  discussion,  the  designations  for  the  real  and 
imaginary  terms  have  been  changed  to  subscript  form.)  The  phase 
velocity,  c,  and  the  attenuation  coefficient,  a,  can  be  obtained  from 
the  complex  velocity: 


1 

c 


Z2 


Re  [ 


Vr+1Vi 


]; 


a 


to 


Im  [ 


J. _ 

VR+ivx 


(3.14) 


In  equation  (3.3),  it  was  established  for  progressive  waves  that 


1_  _  2ac 

Q  a) 


Using  (3.14)  and  (3.15), 

Im  -  =  Re  -  /  2  Q 
v  v  '  x 


(3.15) 


(3.16) 


Rationalizing  the  reciprocal  of  the  complex  velocity  gives 


I  =  VR 
v  vR2+vI2 


VI 

1  — S - 7  ‘ 

vr2+vi2 


From  (3.16)  and  (3.17), 

_  VI 
Vr2+Vi2 

which  gives 


J_  VR 

2Q  (vR2  +  vi2) 


^R 

2<Q 


(3.17) 


vi 


(3.18) 


r.  't.  (  .1  i 


Thus  the  imaginary  part  of  the  complex  velocity  is  related  to  the  real 
part  and  the  specific  attenuation  factor. 
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Although  equation  (3.14)  shows  that  Re  (1/v)  =  1/c,  (3.17)  indi¬ 
cates  that  vr  cannot  be  evaluated  without  prior  knowledge  of  vj,  a 
conflict  with  equation  (3.18).  To  solve  this  problem,  it  can  be  shown 
that  vr  -  c.  From  (3.14),  (3.15)  and  (3.17) 


and 


From  (3. 19) , 


and  solving  for  Vj, 


Re  -  = 


v 


R 


Vj  = 


V 


R 


v  VR2  +  Vj2 


-V- 


Im  -  = 


V  VR2  +  VI2 


1 

c 

1 

2cQ 


VR2  (  VI2 
c  c 


(cvR  -  vr2)1/2 


(3.19) 


(3.20) 


Substituting  into  (3.20), 

-  C  C  Vr  -  Vr2  )l/2 


VR2  +  c  VR  -  Vr2 


1 

2cQ 


Multiplying  by  c  Vr  and  squaring  yields 


v 


C  -  Vr  = 


R 


4  Q< 


Solving  for  vr. 


vR  =  c 


4Q< 


(4Q2  +  1)  • 


Since  Q  >>  1,  the  approximate  relation 

Vr  -  c  (3.21) 

is  established. 

Equation  (3.13)  has  been  programmed  to  generate  the  reflection 
response  of  the  layered  medium  as  a  function  of  frequency.  Some  form  of 
pulse  representing  the  input  waveform  must  also  be  produced.  In  the 
program  written,  the  pulse  normally  used  was  generated  by  setting  its 
amplitude  response  in  the  frequency  domain  to  be  constant  over  a  specified 
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range  of  frequencies.  In  the  time  domain,  this  produces  a  spiked  pulse. 
The  reflection  response  and  Fourier  spectrum  of  the  input  waveform  are 
multiplied  in  the  frequency  domain.  A  fast  Fourier  algorithm  (Gentleman 
and  Sande,  1966)  is  used  to  take  the  inverse  transform  of  the  product, 
this  output  containing  the  amplitude  values  of  the  synthetic  seismogram. 

To  enable  comparison  of  such  a  recording  with  the  observational  data,  a 
facility  was  included  in  the  program  to  obtain  the  autopower  spectrum  of 
any  desired  time  interval  of  the  theoretical  seismogram.  All  outputs 
from  the  various  steps  in  the  program  were  plotted  on  the  CALCOMP  plotter. 
The  Appendix  contains  a  FORTRAN  source  listing  of  the  program  used  to 
generate  synthetic  seismograms  with  spatial  attenuation. 

3.2  Q  of  the  crust 

The  attenuation  of  waves  in  elastic  media  has  been  measured  for 
many  years,  but  it  is  only  within  the  past  half  decade  that  sufficiently 
precise  measurements  and  analytical  techniques  have  been  achieved  to 
enable  computations  of  the  distribution  of  Q  within  the  earth  (Anderson 
and  Archambeau,  1964;  Kovach  and  Anderson,  1964;  Ben-Menahem,  1965  and 
others).  These  studies  generally  rely  on  earthquake-generated  body  and 
surface  wave  data  or  free  oscillations  of  the  earth.  Measurements  of 
attenuation  have  thus  been  made  on  long-period  waves  which  have  travelled 
over  long  distances.  The  resultant  models  of  Q  are  only  reliable  for 
the  mantle  and  not  for  the  crust.  Knopoff  (1964)  and  Sato  (1967)  pro¬ 
vide  comprehensive  reviews  of  the  entire  Q-determination  problem. 

In  contrast  to  the  relatively  large  amount  of  literature  concern¬ 
ing  the  study  of  Q  in  the  upper  mantle,  very  few  results  relating  to  the 
values  of  Q  in  the  crust  have  been  published.  Press  (1964)  determined 
an  average  value  of  Q  =  450  from  the  propagation  of  Lg  waves  from 
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nuclear  explosions  in  Nevada.  His  average  value,  Q  =  260,  for  Pg  waves 
was  lower,  but  this  he  attributed  to  P-wave  energy  loss  by  mode  conver¬ 
sion.  From  a  measurement  of  amplitude  decay  for  Love  and  Rayleigh 
waves  of  periods  50  to  300  seconds,  Anderson  et  al .  (1965)  estimate 
a  Q  for  P-waves  of  about  1000  for  the  crustal  section  to  a  depth  of  38 
kilometers.  At  this  depth,  they  postulate  a  reduction  in  Q  to  a  value 
of  135,  this  Q  extending  for  a  further  22  kilometers.  For  comparison 
with  his  own  observations,  Sumner  (1967)  calculates  a  "grouped”  (or 
average)  Q  of  290  for  Anderson  et  al.’s  model  to  a  depth  of  180  kilo¬ 
meters.  He  points  out  the  disagreement  with  his  Q  of  greater  than 
1000  for  frequencies  from  1.5  to  15  Hz.  and  depths  from  the  surface 
to  180  kilometers.  This  Q  was  established  from  a  study  of  local 
earthquakes  in  southern  Peru.  From  a  study  of  shear  waves  propagating 
in  the  uppermost  mantle  and  for  paths  distributed  throughout  the  world, 
Molnar  and  Oliver  (1969)  show  there  is  a  significant  lateral  variation 
of  upper  mantle  attenuating  properties.  They  demonstrate  that  the 
shear  waves  are  propagated  efficiently  throughout  the  stable  regions 
of  the  earth,  but  inefficiently  through  active  regions  such  as  island 
arcs  and  ocean  ridges.  Thus  reported  values  of  Q  can  be  highly 
dependent  on  the  regions  in  which  the  seismic  measurements  were  made. 

In  the  present  research,  an  attempt  is  made  to  derive  a  Q- 
structure  for  the  continental  crust  in  southern  Alberta  from  an 
analysis  of  the  autopower  spectra  of  reflections  from  the  deep  crust 
and  the  results  from  theoretical  seismograms  incorporating  spatially 
dependent  attenuation.  A  description  of  the  autopower  spectra  for 
deep  reflections  on  a  number  of  seismograms  was  included  in  Section 
2.2.  Such  spectra  were  computed  over  traveltime  intervals  including 
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the  R  and  M  reflections  for  most  available  digitized  seismograms. 

One  characteristic  which  was  noted  for  many  of  these  spectra  was  a  gen¬ 
eral  trend  toward  less  power  at  the  higher  frequencies  within  the  5  to 
25  Hz.  interval  studied.  Assuming  that  the  trend  was  due  primarily  to 
the  intrinsic  attenuation  of  seismic  energy  as  it  propagated  through  the 
crustal  rocks,  it  was  decided  to  examine  the  rate  of  such  attenuation. 

In  order  to  ensure  that  reflected  energy  was  being  considered,  only  those 
autopower  spectra  for  which  the  quality  of  the  reflections  in  the  lower 
crust  at  two  separate  time  intervals  was  fairly  good  were  used  in  the 
analysis.  Figure  3.1  illustrates  several  autopower  spectra  for  a  two- 
second  interval  centered  about  the  R  or  M  reflection  times.  It  should 
be  noted  that  the  instrument  responses  have  been  removed  from  these  curves. 
The  straight  lines  drawn  across  the  individual  spectra  indicate  the  general 
trend  toward  less  power  for  the  higher  frequency  components  of  the  re¬ 
flected  wave  trains.  To  the  left  of  each  curve,  the  slope  of  the  line  is 
given.  The  sole  purpose  in  presenting  this  linear  relationship  is  to  pro¬ 
vide  a  means  of  comparing  the  rate  of  attenuation  on  the  spectra  of  observed 
reflected  wavelets  with  similar  quantities  determined  from  the  autopower 
spectra  of  synthetic  seismograms.  This  method  of  analysis  presupposes  a 
flat  source  spectrum  over  the  frequency  interval  studied.  The  examples 
were  chosen  to  show  the  manner  in  which  the  slopes  are  measured  and  to 
illustrate  typical  variations  in  the  autopower  spectra.  Record  39,  1966 
illustrates  one  of  the  lower  slope  values  measured;  Record  32,  1967  pro¬ 
vides  an  example  of  one  of  the  higher  calculated  rates.  The  attenuation 
rates  of  Records  01  and  45,  1966  represent  more  usual  values.  By  making 
such  measurements  on  many  autopower  spectra,  an  indication  of  the  char¬ 
acteristic  values  is  obtained.  Figure  3.2  is  a  histogram  of  the  rates 
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Figure  3.1 


Autopower  spectra  for  selected  time  intervals 
on  four  field  seismograms.  The  straight  lines 
indicate  the  rates  at  which  the  frequency 
components  are  being  attenuated. 
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of  attenuation  for  data  from  52  seismograms.  The  figure  shows  that  most 
of  the  measurements  give  values  ranging  from  0.30  to  0.55  db/Hz.  as  the 
basic  rate  at  which  power  is  being  lost  in  the  higher  frequencies.  At 
the  bottom  of  the  diagram,  the  average  slope  for  the  M  and  R  traveltime 
intervals  is  given.  The  two  rates  of  attenuation  are  nearly  identical, 
being  0.42  and  0.41  db/Hz.  for  the  M  and  R  reflections,  respectively. 

The  method  described  in  Section  3.1  was  used  to  generate  the 
reflection  coefficients  as  a  function  of  frequency  for  a  horizontally 
layered  earth  with  depth-dependent  attenuation.  A  velocity  log,  for 
which  the  values  were  averaged  over  20  foot  intervals,  from  a  well 
(Tenneco  Eyremore  10-15,  Lsd.10,  Sec. 15,  Twp.18,  Rge.19  W4M)  in  the 
region  of  the  reflection  profiles  was  used  to  model  the  sedimentary 
strata  to  a  depth  of  1.68  kilometers,  the  extent  of  the  log.  Because 
the  method  used  does  not  generate  reflections  at  the  free  surface,  a  near¬ 
surface  velocity  contrast  of  1.22  to  2.29  km/sec  (4000  to  7500  ft/sec) 
was  introduced  to  produce  the  multiple  reflections  from  this  boundary. 

The  remainder  of  the  crustal  section  was  defined  by  a  transition  zone 
at  the  depths  of  the  R  (34  kilometers)  and  M  (47  kilometers)  discontin¬ 
uities.  In  the  first  case,  this  was  considered  to  be  a  first-order 
discontinuity  in  which  the  velocities  abruptly  changed  from  6.5  to  7.3 
km/sec  for  R  and  from  7.4  to  8.3  km/sec  for  M,  these  values  being  inferred 
from  refraction  data.  Specification  of  the  Q  structure  for  the  desired 
crustal  model  completed  the  input  data  necessary  to  determine  the 
crustal  reflection  response. 

It  is  also  necessary  to  specify  an  input  waveform  in  order  to 
generate  the  synthetic  seismogram  for  plane  waves  at  normal  incidence. 
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Figure  3.2 


Histogram  of  the  rates  of  attenuation  of 
frequency  as  measured  from  52  seismograms . 


Number  of  seismograms 


12 


i  I 
I  I 
I 
I 


- 

1 

1 

1 

1_ L 

1 

l 

l 

i 

! 

_ 

0 

o 

CNJ 

.40 

.60 

• 

00 

o 

Rate  of  attenuation  of  frequency 
(  db  per  Hz.  ) 


M  (52  seismograms):  average  =  0.42  db/Hz. 
R  (52  seismograms):  average  =  0.41  db/Hz. 


:  (  r  v  r 


70 


The  frequency  content  of  energy  generated  by  an  explosion  in  a  borehole 
is  not  well  known.  Some  indication  of  the  spectral  characteristics  from 
an  explosive  source  in  a  homogeneous  medium  may  be  obtained  from  studies 
of  Sharpe  (1942a  and  1942b)  and  Heelan  (1953).  However,  the  lower  fre¬ 
quencies  have  probably  been  enhanced  over  the  theoretical  case  in  our 
field  experiments  by  the  presence  of  inhomogeneities  and  also  the  use 
of  a  pattern  of  five  holes.  From  a  study  of  observations  made  near  a 
small  (1/8  lb.)  charge,  Sharpe  (1942b)  notes  that  even  in  competent 
material,  there  is  an  observable  dissipation  of  frequencies  as  low  as 
150  Hz.  at  distances  of  only  60  meters.  In  petroleum  exploration,  the 
frequencies  of  the  reflected  energy  are  usually  in  a  range  from  30  to 
80  Hz.  (Jakosky  and  Jakosky,  1952).  The  deep  reflection  data  from  this 
study  and  others  is  comprised  mainly  of  frequencies  from  7  to  30  Hz. 
Since  we  had  no  prior  knowledge  of  the  relative  amplitudes  of  different 
frequency  components  generated  by  the  explosions,  it  was  decided  to 
produce  a  pulse  input  which  had  a  constant  amplitude  spectrum  from  7 
to  70  Hz.  The  remainder  of  the  spectrum  was  specified  by  a  decrease 
to  zero  amplitude  from  7  to  0  Hz.  and  70  to  77  Hz.  Figure  3.3  shows 
the  input  pulse  resulting  from  such  frequency  specifications  together 
with  its  amplitude  spectrum. 

The  general  character  of  the  transfer  function  of  the  layered 
media  is  included  in  Figure  3.3,  although  the  details  change  depending 
upon  the  particular  model  chosen.  An  inverse  Fourier  transform  of  the 
product  of  the  complex  spectra  of  the  incident  pulse  and  layered  medium 
produces  the  synthetic  seismogram  shown  on  the  lower  right  part  of  the 
figure.  A  variable  gain  control  factor  has  been  applied  only  to  the 
plotted  seismogram  to  provide  reasonable  amplitude  values  for  display. 
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Figure  3.3  Various  quantities  computed  by  the  synthetic 

seismogram  program .  An  inverse  Fourier  trans¬ 
form  of  the  product  of  the  complex  spectrum 
of  the  input  pulse  and  the  transfer  function 
of  the  layered  media  generates  the  synthetic 
seismogram .  The  autopower  spectrum  of  the 
latter  for  the  interval  0.4  to  3.4  seconds  is 
also  given . 
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In  the  upper  right  section  of  the  illustration,  the  autopower  spectral 
density  of  the  synthetic  seismogram  from  0.4  to  3.4  seconds  is  plotted. 
For  all  calculations  of  autopower  spectra,  the  amplitudes  not  modified 
by  the  gain  control  were  used.  Since  only  the  Q  values  of  the  upper 
section  were  varied  in  the  following  investigations,  this  spectrum  is 
representative  of  the  power  estimates  of  reflections  from  the  sediments 
for  any  case  considered. 

The  first  model  chosen  consisted  of  first-order  discontinuities 
at  depths  of  the  R  and  M  interfaces.  Different  Q  structures  were 
attempted  to  determine  the  effects  of  this  attenuation  on  the  autopower 
spectrum  of  the  reflections  on  the  synthetic  seismograms.  For  an 
effective  Q  of  infinity  throughout  the  crust,  that  is  no  attenuation, 
the  general  trend  showed  no  decrease  of  energy  content  with  higher 
frequencies.  Finite  values  of  Q  were  then  introduced  and  it  was  soon 
discovered  that  the  dominant  trend  for  frequency  attenuation  was  con¬ 
trolled  by  Q  in  the  sedimentary  layers. 

Figure  3.4  illustrates  the  autopower  spectra  for  reflections 
from  the  R  and  M  horizons  modelled  as  first-order  discontinuities  for 
three  different  Q  structures  in  the  sediments  with  an  infinite  Q  below. 
The  depth  extent  to  which  each  Q  applied  was  set  to  correspond  with 
the  major  velocity  changes  observed  on  the  continuous  velocity  log. 

The  first  Q,  extending  over  only  5  meters,  was  included  to  model  the 
effects  of  the  surface  layer.  A  comparison  of  the  rates  at  which 
higher  frequency  components  are  being  attenuated  can  be  made  between 
the  histogram  of  Figure  3.2  and  the  slopes  measured  and  noted  in 
Figure  3.4.  The  distributions  of  Q  represented  by  the  models  in  the 
upper  two  figures  generate  rates  of  attenuation  which  are  lower  than 
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Figure  3.4 


Effects  of  different  distributions  of  Q  on 
the  rates  of  attenuation  of  frequency  for 
reflections  centered  about  the  R  and  M 
horizons  modelled  as  first  order  discontinuities . 
R  discontinuity:  6.50  to  7.31  km/sec  at  a 

depth  of  34  kilometers . 

M  discontinuity:  7.36  to  8.25  km/sec  at  a 

depth  of  47  kilometers. 

The  variation  of  Q  with  depth  for  each  model 
is  shown  in  the  upper  right  comer  of  each 
plot.  The  rates  of  attenuation  in  db/Hz .  are 
also  indicated. 
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those  generally  observed.  For  the  model  of  the  lower  left  part  of  the 
figure,  for  which  Q  values  are  one  half  those  in  the  upper  right,  it  is 
seen  that  the  slopes  have  become  higher  than  the  average  given  on  the 
histogram.  In  the  lower  right  diagram,  an  upper  Q  structure  identical 
with  the  illustration  above  it  has  been  assumed  but  a  constant  value 
of  Q  =  1500  has  been  included  from  the  base  of  the  sediments  to  the  top 
of  the  mantle.  This  causes  an  increase  in  the  rate  of  attenuation  for 
both  R  and  M  reflection  groups.  A  relative  change  in  the  rates  between 
R  and  M,  the  values  of  which  were  nearly  identical  for  no  attenuation  in 
the  lower  crust,  is  also  noted.  But  the  histogram  shows  evidence  indi¬ 
cating  that  the  attenuation  rates  are  about  the  same.  This  would 
suggest  that  the  value  of  Q  between  the  R  and  M  groups  of  reflections 
should  be  increased.  From  the  preceding  discussion,  a  possible  Q 
structure  has  an  average  value  of  290  in  the  sedimentary  rocks,  the 
crystalline  basement  has  a  Q  of  about  1500  and  the  lower  crust  would 
appear  to  have  a  Q  in  excess  of  1500. 

The  only  example  of  Q  values  for  a  shallow  section  is  that  given 
in  a  hypothetical  case  by  Trorey  (1962)  to  exemplify  the  effects  of 
absorption  on  the  first  2.5  seconds  of  exploration  seismograms.  The 
values  chosen  in  this  study  are  not  inconsistent  with  Trorey* s  if  one  takes 
into  account  the  generally  higher  velocities  and  thus  higher  Q  values  in 
the  southern  Alberta  section.  A  comparison  in  the  time  domain  of  the 
synthetic  seismogram  from  the  sedimentary  layers  with  a  field  record 
from  the  area,  obtained  through  the  courtesy  of  Tenneco  Oil  and  Minerals, 
Ltd. ,  indicates  a  reasonable  degree  of  qualitative  agreement. 

One  point  which  follows  from  the  preceding  discussion  is  the 
necessity  of  having  wide-band  recording  of  reflections  from  the  post- 
Precambrian  layers.  If  power  spectra  of  such  data  could  be  obtained, 
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then  more  precise  Q  models  for  the  sedimentary  strata  could  be  achieved. 
This  would  enable  Q's  for  the  lower  crust  to  be  specified  with  more 
reliability.  It  must  also  be  recalled  that  the  foregoing  results  were 
based  on  an  assumption  of  a  single  first-order  discontinuity  for  both 
the  R  and  M  transition  zones.  It  is  shown  in  the  next  section  that 
such  a  model  does  not  agree  with  observational  results.  When  more  complex 
velocity  models  for  the  transition  zones  are  considered,  it  is  found  that 
by  their  very  nature  such  zones  may  selectively  enhance  or  decrease 
certain  frequency  components,  thus  introducing  an  apparent  attenuation. 

3.3  Properties  of  the  crustal  reflectors 

On  the  basis  of  the  results  in  Germany  (Dohr  and  Fuchs,  1967),  in 
the  Soviet  Union  (Beloussov  et  al.,  1962)  and  in  the  present  study,  it 
is  well  established  that  the  deep  crust  is  reflecting  seismic  energy.  A 
question  concerning  the  properties  of  these  reflecting  zones  within  the 
crustal  section  can  be  justifiably  raised.  Fuchs  (1969)  has  recently 
pointed  out  that  the  classical  model  of  a  simple  layered  crust  is  not 
concordant  with  a  number  of  observations  from  the  deep  reflection  work 
in  Germany.  Such  observations  include  reflections  which  rarely  correlate 
over  more  than  a  few  kilometers,  the  band  character  of  the  reflected 
signals,  the  apparent  existence  of  a  lower  cut-off  frequency  and  the 
recording  of  amplitudes  larger  than  would  be  expected  from  a  first-order 
discontinuity.  To  explain  these  observations,  Fuchs  postulates  that 
the  reflecting  regions  are  best  represented  by  laminated  transition 
zones  in  which  there  is  a  series  of  abrupt  velocity  increases  and 
decreases  over  small  depth  intervals  of  less  than  120  meters.  He  shows 
that  the  reflection  coefficients  for  such  models  are  relatively  large 
over  certain  ranges  of  frequency  and  that  the  reflected  wavelets  do 
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exhibit  a  band  character. 

A  few  remarks  can  be  made  concerning  the  deep  reflection  observa¬ 
tions  which  Fuchs  used  to  postulate  such  transition  zones.  We  have 
previously  shown  that  some  deep  crustal  reflections  may  be  continuously 
correlated  over  many  tens  of  kilometers  (Clowes  et  al.,  1968).  In  the 
previous  chapter,  a  comment  was  made  concerning  the  generation  of 
seismograms  with  an  oscillatory  appearance  when  the  original  data  were 
filtered  by  a  narrow  band-pass  filter.  Finally,  it  might  be  noted  that 
part  of  the  reason  that  near-vertical  incidence  reflections  with  fre¬ 
quencies  below  5  Hz.  have  not  been  successfully  recorded  might  be  due 
to  source  conditions.  Sharpe  (1942a)  notes  that  an  important  frequency 
component  generated  from  the  application  of  an  explosive  pressure 
within  a  spherical  cavity  in  an  elastic  medium  will  be 

f  -  JiL  -  —  Z__ 

2tt  3tt  Rs 

where  v  is  the  velocity  of  the  medium  and  Rs  is  the  radius  of  the  cavity. 
This  implies  that  to  produce  a  dominant  component  of  5  Hz.  in  a  medium 
of  velocity  1  km/sec,  the  cavity  radius,  which  in  a  real  medium  may  be 
approximately  equivalent  to  the  radius  at  which  the  stress  produced 
equals  the  strength  of  the  material,  would  be  about  30  meters.  To  gener¬ 
ate  lower  frequencies  would  require  a  medium  with  unusually  low  velocities 
or  a  very  large  effective  cavity.  Thus  it  is  possible  that  the  small 
explosions  typically  used  in  boreholes  may  not  generate  much  energy  at 
the  lower  frequencies. 

The  synthesis  of  seismograms  with  attenuation  as  described  in 
Section  3.1  has  been  used  to  provide  a  means  by  which  transition  zones 
can  be  studied.  Autopower  spectra  of  selected  time  intervals  on  such 
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data  may  be  computed  and  compared  with  those  calculated  from  the  field 
seismograms.  Figure  3.5  illustrates  the  autopower  spectral  densities 
of  theoretical  seismograms  over  a  three-second  interval  for  three  dif¬ 
ferent  types  of  transition  zones  at  the  depths  of  the  R  and  M  reflecting 
horizons.  For  display  purposes,  the  spectrum  for  M  in  models  B  and  C 
has  been  raised  by  about  one  order  of  magnitude  above  that  for  R.  The 
characteristics  of  the  assumed  Q  distribution,  the  type  of  transition 
zones  and  the  form  of  the  reflected  wavelet  for  each  model  are  given 
in  Figure  3.6.  Note  that  the  peak  values  of  the  spectra  for  models 
B  and  C  are  approximately  one  and  two  orders  of  magnitude,  respectively, 
greater  than  the  maximum  value  for  model  A.  This  has  an  important 
bearing  on  the  acceptance  of  a  proposed  model  for  it  relates  to  the 
amplitudes  of  the  reflected  wavelets.  From  the  synthetic  seismogram 
for  the  sedimentary  layers,  the  average  peak  power  in  the  interval  0.4 
to  3.4  seconds  is  about  102-7  in  relative  units  (Figure  3.3).  It  should 
be  recalled  that  the  synthetic  seismograms  generated  are  for  plane  waves 
at  normal  incidence  in  a  horizontally  layered  medium.  Thus  the  spheri¬ 
cal  spreading  factor  (1/r)  is  not  taken  into  account.  The  dominant 
reflector  within  the  sedimentary  section  has  a  two-way  depth  of  three 
kilometers  while  the  two-way  depths  to  the  R  and  M  interfaces  are 
approximately  60  and  80  kilometers.  Assuming  a  decrease  of  power  at  a 
rate  of  1/r2,  it  is  found  that  the  R  and  M  reflections  are  attenuated 
by  factors  of  10  2,6  and  lO"2*9,  respectively,  relative  to  the  shallow 
reflector. 

During  the  course  of  the  field  program,  a  few  attempts  were  made 
to  obtain  reflections  from  the  sedimentary  layers  with  the  same  record¬ 
ing  procedure  used  for  the  deep  reflection  work,  except  that  only  one 
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Figure  3„5 


Autopower  spectra  for  three-second  time 
intervals  centered  about  the  R  and  M 
reflections  on  synthetic  seismograms „  The 
ordinate  scales  are  in  arbitrary  units „ 

In  ($)  and  (C)3  the  curve  for  M  has  been 
raised  by  about  one  order  of  magnitude  for 
display  purposes o  The  types  of  model  for 
the  R  and  M  horizons  and  the  Q- structures 
which  were  used  in  generating  the  theoretical 
seismograms  are  given  in  Figure  30  6e 
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Figure  3.6 


Details  of  the  models  for  which  the  auto¬ 
power  spectra  of  Figure  3„5  are  computedo 
The  shape  of  the  wavelet  reflected  from 
the  indicated  transition  zone3  for  a  pulse 
input  waveform  (Fig-ure  3.3)  3  is  also 
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five-pound  charge  in  one  hole  was  used.  These  attempts  were  generally 
unsuccessful  because  the  field  techniques  were  not  optimized.  However, 
in  one  instance,  a  seismogram  with  identifiable  reflections  from  about 
the  top  of  the  Paleozoic  and  one  deeper  horizon  was  recorded.  A  normal 
deep  reflection  seismogram  was  also  recorded  for  the  same  positions 
of  the  detectors.  A  comparison  of  the  relative  amplitudes  of  the 
shallow  and  deep  reflection  events  could  thus  be  made.  When  the  data 
were  corrected  for  different  amplifier  attenuation  settings  and  assuming 
a  linear  dependence  of  resultant  amplitude  on  charge  size,  a  ratio  of 
53:1  was  calculated.  Squared  to  provide  an  estimate  of  relative  power, 
this  ratio  becomes  approximately  103o4:l.  Such  a  figure  can  be  used  as 
a  guide  in  evaluating  the  applicability  of  the  ratio  of  power  spectra 
for  the  shallow  and  deep  reflections  generated  in  the  synthetic  seismo¬ 
gram,  after  allowance  has  been  made  for  spherical  spreading. 

From  the  characteristics  of  the  observed  autopower  spectra  of 
Figure  3.1  and  the  spectra  for  the  case  of  first-order  discontinuities 
in  Figure  3.4,  it  is  evident  that  the  latter  form  of  transition  zone 
does  not  yield  spectral  characteristics  similar  in  form  to  those 
observed.  Thus  more  complex  models  had  to  be  constructed.  Model  A  of 
Figure  3.5  illustrates  the  spectral  densities  for  a  zone  with  a  linear 
increase  of  velocity  over  0.5  kilometers.  This  produces  the  two  peaks 
often  observed  but  even  with  the  model  incorporating  a  distribution  of 
relatively  high  Q's,  the  rates  of  attenuation  of  the  higher  frequency 
components  (.70  and  1.0  db/Hz.  for  R  and  M,  respectively)  are  much 
larger  than  the  observed  values .  The  shape  of  the  reflected  wavelet 
for  a  pulse  input  bears  no  relation  to  any  observed  reflections. 

Because  the  peak  power  of  the  curve  for  R  is  about  10-0°6,  the  ratio 
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of  energy  in  the  shallow  reflections  to  that  in  the  R  reflections  is 
105‘9:1,  after  allowance  has  been  made  for  the  effects  of  spherical 
spreading.  This  compares  to  the  observed  ratio  of  10 3  *  4 : 1 .  It  can  be 
noted  that  a  gradient  over  only  0.25  kilometers  was  also  tried,  but  the 
shape  of  the  curve,  rate  of  attenuation  and  relative  amplitudes  did  not 
agree  with  observation.  For  these  reasons,  the  possibility  of  a  linear 
increase  in  velocity  with  depth  as  a  model  for  transition  zones  can 
probably  be  eliminated. 

The  second  transition  model,  B,  of  Figure  3.5  comprises  three  equal 
step  increases  in  velocity,  the  total  extent  of  the  zone  being  0.8  kilo¬ 
meters.  The  general  shapes  of  the  computed  spectra  correlate  well  with 
those  for  some  observed  data.  Although  the  measured  rates  of  attenuation 
(.60  and  .77  db/Hz.  for  R  and  M,  respectively)  are  somewhat  high,  the 

i 

qssumed  distribution  of  Q  was  low  so  that  attenuation  was  considerable. 

The  pulse  input  reflected  from  such  a  zone  produces  a  sequence  of  pulses 
with  strong  amplitudes  in  one  direction  only.  However,  this  wavelet  form 
bears  a  closer  resemblance  to  observed  reflections  than  does  that  of 
model  A.  Comparing  the  power  in  the  shallow  reflections  with  the  peak 
value  for  the  R  curve,  and  considering  the  effects  of  spherical  spreading, 
the  ratio  is  104°8:1  which  is  more  than  one  order  of  magnitude  greater 

than  that  measured  from  the  recorded  seismograms.  If  the  depth  extent  of 
the  transition  zone  is  lessened,  everything  else  remaining  constant,  the 
effect  is  a  shift  of  the  frequencies  at  which  the  spectral  peaks  occur, 
but  the  slopes  and  power  spectrum  ratio  are  relatively  unaffected.  Thus 
a  step  transition  zone  as  described  by  model  B  can  only  be  considered  to 
have  a  low  probability  of  modelling  the  deep  reflectors. 

Model  C  is  based  on  the  preferred  solution  for  the  deep  reflecting 
zones  as  described  by  Fuchs  (1969).  It  consists  of  a  series  of  sharp 
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increases  and  decreases  of  velocity,  for  which  the  thickness  of  the 
individual  layers  is  180  meters.  The  shapes  of  the  spectral  curves 
bear  a  good  similarity  to  many  of  those  computed  from  the  deep  reflec¬ 
tion  seimograms.  As  measured  from  the  spectrum,  the  rates  of  attenua¬ 
tion  of  frequency  for  the  R  and  M  curves  are  .52  and  .53  db/Hz., 
respectively.  These  rates  are  in  the  range  of  those  given  in  the 
histogram  of  Figure  3.2.  It  is  seen  that  the  attenuation  included 
in  this  case  is  relatively  high  since  the  Q’s  are  low.  From  the 
shape  of  the  reflected  wavelet,  it  is  evident  that  the  transition 
region  of  model  C  produces  an  oscillatory  sequence  of  reflected  pulses 
from  a  single  input  pulse.  The  ratio  of  the  power  in  the  spectrum  of 
the  shallow  reflections  to  the  peak  power  in  that  of  the  R  reflection, 
allowing  for  the  effects  of  spherical  spreading,  is  103*5:1.  This 
agrees  favorably  with  the  calculated  ratio  of  1 0 3  * 4 : 1  for  observed 
reflections  on  the  field  seismograms.  The  model  of  a  reflecting 
interface  consisting  of  thin,  multiple  layers  of  high  and  low  velo¬ 
cities  best  satisfies  all  the  observational  data. 

In  Figure  3.7,  a  comparison  of  the  autopower  spectral  densities 
for  a  two-second  interval  centered  about  the  R  and  M  reflection  times 
of  a  recorded  seismogram  (solid  curves)  is  made  with  the  spectra 
computed  from  a  synthetic  seismogram  (dashed  curves).  The  transition 
model  for  both  the  R  and  M  reflecting  zones  and  the  specified  structure 
of  Q  are  given  in  the  lower  part  of  the  figure.  The  agreement  between 
the  observational  and  theoretical  curves  is  quite  good.  In  terms  of 
the  frequencies  at  which  the  spectra  have  their  peak  values,  the  depth 
extent  of  the  individual  high  and  low  velocity  zones  is  critical.  For 
example,  if  the  depth  interval  of  the  R  zone  was  reduced  to  0.12  from 
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Figure  3.7  Comparison  of  the  autopower  spectra  for 

two-second  intervals  centered  about  the 
R  and  M  reflections  on  field  seismograms 
(solid  curves)  with  spectra  computed  from 
reflections  on  synthetic  seismograms 
(dashed  curves ),  The  ordinate  scale  is 
in  relative  units .  The  variation  of  Q 
with  depth  and  the  characteristics  of  the 
transition  zone  models  are  shown  in  the 
lower  part  of  the  figure. 
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0.15  kilometers,  the  resultant  spectrum  would  peak  at  about  15  Hz. 

A  variation  in  the  individual  thicknesses  of  the  velocity  lenses  can 
thus  account  for  the  observed  variations  in  the  peak  frequencies  of 
the  spectra.  However,  it  should  be  emphasized  that  the  total  extent 
of  the  transition  zone  is  less  than  one  kilometer. 

A  slightly  more  complex  version  of  the  laminated  model,  having 
thicknesses  of  0.1  and  0.2  kilometers  for  the  high  and  low  velocity 
layers,  respectively,  has  been  attempted.  This  yielded  a  spectrum 
for  which  the  frequency  of  the  maximum  power  was  similar  to  that  for 
equal  layers  of  0. 15-kilometer  thickness.  The  width  of  the  spectral 
peak  was  somewhat  broadened  relative  to  the  case  of  equal  layers. 

Similar  models  could  probably  be  designed  to  produce  spectra  which 
would  closely  approach  observed  spectral  densities  with  broad  maxima. 
Another  version,  in  which  gradients  over  0.2  kilometers  replaced  the 
abrupt  velocity  increases  and  decreases,  resulted  in  a  spectrum  which 
had  slightly  lower  amplitudes,  a  high  rate  of  attenuation  and  a  form 
not  very  similar  to  the  observed  ones.  One  further  effect  was  investi¬ 
gated.  Two  sets  of  identical  lenticular  transition  zones  were  modelled 
to  produce  two  sequences  of  reflected  pulses  separated  by  about  one 
second  of  two-way  traveltime.  The  first-order  multiples  generated 
between  these  zones  had  negligible  amplitudes  on  the  synthetic  seismo¬ 
gram.  An  autopower  spectrum  over  both  reflected  wavelets  was  calculated 
and  showed  that  all  the  basic  features  of  the  spectrum  from  just  one  of 
the  pulse  sequences  was  retained.  Irrespective  of  the  particular  details 
for  any  one  model,  it  was  found  that  some  form  of  transition  zone  which 
included  velocity  reversals  over  small  depth  intervals  was  necessary 
to  produce  results  concordant  with  observations. 
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Because  the  autopower  spectra  of  the  two-second  time  intervals, 
centered  about  reflections  interpreted  as  being  from  the  M  discontinuity, 
include  effects  of  crustal  properties  over  a  depth  extent  of  about  seven 
kilometers,  these  spectra  do  not  represent  just  the  M  itself.  In 
Chapter  4,  it  is  shown  that  a  number  of  reflections  appear  to  come  from 
within  the  deepest  part  of  the  crust,  so  that  much  of  the  character  of 
the  computed  spectra  from  the  field  seismograms  might  be  attributed  to 
effects  near  the  base  of  the  crust.  Thus  the  form  of  the  M  reflector 
illustrated  in  Figure  3.7  does  not  necessarily  represent  the  Mohorovicid 
discontinuity,  but  may  model  a  number  of  reflecting  zones  at  depths 
near  the  Moho.  However,  it  can  be  reiterated  that  on  the  basis  of  a 
detailed  study  of  wide  angle  reflections  from  the  M  discontinuity, 
Meissner  (1967a)  proposed  a  laminated  structure  with  individual  layers 
less  than  0.15  kilometers  thick. 

The  model  of  a  deep  crustal  reflector  which  has  alternating  lenses 
of  low  and  high  velocity  material  poses  the  problem  of  how  such  a  zone 
could  be  formed  by  natural  processes.  It  could  possibly  represent  a 
series  of  intermittently  occurring  lava  flows  in  which  sediments,  since 
metamorphosed,  were  deposited  during  periods  of  quiescence.  Perhaps  it 
represents  the  effects  of  partial  melting  and  the  segregation  of  more 
basic,  high  velocity  material  from  lower  velocity,  acidic  rocks.  From 
petrological  studies,  Ringwood  and  Green  (1966)  give  evidence  that 
mineral  assemblages,  which  are  stable  under  the  pressure-temperature 
conditions  of  the  lower  continental  crust  and  which  are  likely  to  exist 
there,  have  densities  and  seismic  velocities  that  are  much  higher  than 
values  believed  to  be  characteristic  of  the  lower  crust.  They  suggest 
that  in  this  region  there  exist  large  quantities  of  minerals  of 
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relatively  low  seismic  velocity  to  counterbalance  the  higher  velocities 
of  the  other  rocks.  Such  a  hypothesis  could  have  some  relation  to  the 
postulated  deep  crustal  reflectors. 
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CHAPTER  4 

INTERPRETATION  OF  THE  DATA 
4.1  The  analog  reflection  seismograms 

In  the  last  chapter,  it  was  demonstrated  that  inherent  on  the 
recorded  seismograms  was  information  relating  to  the  nature  of  the  reflec¬ 
ting  horizons  and  to  the  elastic  properties  of  the  crustal  section. 
However,  the  primary  intent  of  the  experimental  research,  and  the  philo¬ 
sophy  behind  the  design  of  the  project,  was  an  investigation  of  possible 
structural  features  in  the  lower  crust.  To  this  end,  more  than  one 
hundred  deep  reflection  seismograms  from  the  Lomond  and  Blackfoot  profiles 
were  carefully  processed  and  analysed  to  provide  as  comprehensive  an 
interpretation  as  possible.  Since  a  considerable  time  lag  often  existed 
between  the  data  acquisition  in  the  field  and  the  subsequent  digitization 
in  the  laboratory,  an  analysis  of  the  photographically  recorded  seismo¬ 
grams  was  considered  necessary  and  important.  Usually  the  quality  of 
these  data  was  sufficiently  good  to  enable  a  preliminary  interpretation 
of  the  lower  crustal  structure. 

For  such  an  interpretation,  the  seismograms  were  examined  to 
identify  most  coherent  phases  which  could  possibly  be  identified  with 
reflected  energy.  Because  the  field  method  employed  was  continuous  pro¬ 
filing,  two  traces  on  adjacent  records  had  reflected  energy  which  ideally 
followed  the  same  travel  path  but  in  opposite  directions,  thus  allowing 
a  time  correlation  (time-tying)  of  reflected  events  between  these  records. 
Phase  correlation  of  the  seismograms  also  proved  possible.  In  this 
manner,  many  reflecting  events  were  correlated,  some  over  much  of  the 
entire  profile,  others  over  only  a  few  kilometers.  A  preliminary  time 
cross  section  of  the  profile  was  plotted,  with  each  reflection  being 
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placed  directly  below  the  surface  region  at  which  it  was  recorded.  From 
the  steep  dips  which  could  be  seen,  it  was  evident  that  all  reflections 
would  have  to  be  migrated  in  order  to  display  the  reflecting  horizons  in 
their  true  subsurface  positions.  For  this  purpose,  a  wavefront  migration 
chart  was  generated  following  the  method  described  by  Gaby  (1945).  A 
velocity  function  consisting  of  five  layers  of  constant  velocity  was  used. 
In  the  post- Pre Cambrian  layers,  the  velocities  were  determined  from  con¬ 
tinuous  velocity  logs  of  wells  in  the  vicinity,  while  refraction  and 
reflection  data  were  utilized  in  the  lower  section  (Cumming  and  Kanasewich, 
1966).  In  Figure  4.1,  the  migrated  seismic  cross  section  is  plotted  as 
a  function  of  two-way  traveltime  on  an  approximately  1:1  scale.  A  depth 
scale  has  been  appended  to  the  section  to  aid  interpretation  but  it  must 
be  remembered  that  along  the  section  the  depth  scale  varies  due  to  varia¬ 
tions  in  average  velocity  as  the  thickness  of  the  layers  changes.  The 
depth  on  the  right  applies  to  the  deepest  part  of  the  structure;  that 
near  the  middle  applies  to  the  highest  part. 

The  lack  of  reflections  between  0  and  6  to  9  seconds  is  due  to  a 
combination  of  factors.  As  mentioned  previously,  the  instrumentation 
incorporates  no  automatic  gain  control  so  that  the  signals  received  between 
0  and  4  seconds  generally  exceeded  the  dynamic  range  of  the  amplifiers. 
Additionally,  large  amplitude  surface  waves  generated  by  the  explosion 
often  obscure  reflected  energy  up  to  times  of  6  or  9  seconds.  As  evi¬ 
denced  by  the  figure,  this  problem  was  more  severe  for  the  southern  pro¬ 
file  than  for  the  northern  profile,  where  coherent  events  from  shallower 
events  were  observed. 

The  reflection  labelled  R  was  the  one  with  the  largest  amplitude 
on  most  seismograms  and  appeared  to  correlate  across  nearly  the  entire 
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Figure  4.1  (a)  and  (b) :  Residual  total  magnetic  intensity 

and  Bouguer  gravity respectively 3  along  the 
seismic  profiles .  (c)  Surface  topography  along 

the  profiles .  The  vertical  exaggeration  is 
about  ten  times .  Arrows  indicate  shot  point 
locations .  (d)  Reflection  seismic  cross  section , 

depths  and  velocities  of  refracting  horizons  (**) 
from  an  east-west  survey  are  superimposed  at  the 
position  of  intersection .  T/ze  "*6 .  2  km/sec " 
velocity  is  an  average  vertical  velocity  between 
the  top  of  the  Paleozoic  and  the  Riel  (R) 
continuity  as  obtained  from  reflection  data . 
Correlation  of  the  reflecting  events  is  indicated 
by  dots . 

Rote  the  east-west  offset  of  16  kilometers  between 
the  two  halves  of  the  section „ 
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90-kilometer  profile .  At  the  position  where  a  reversed  refraction  survey 
had  been  made  (Camming  and  Kanasewich,  loc.  cit.)  the  depth  to  this  re¬ 
flector  corresponded  with  the  depth  of  an  intermediate  refracting  horizon. 
Previously  this  interface  was  termed  the  "Conrad"  discontinuity  (Richards 
and  Walker,  1959;  Kanasewich  and  Cumming,  1965).  This  nomenclature  has 
been  criticized  as  it  could  imply  a  stratigraphic  or  petrologic  correla¬ 
tion  with  the  European  Conrad  discontinuity.  The  latter  represents  an 
interface  with  velocity  changing  from  about  5.6  to  6.2  km/sec  (Conrad, 
1925;  Jeffreys,  1926).  In  southern  Alberta,  the  refracting  and  reflecting 
intermediate  layer  represents  a  change  in  velocity  from  6.5  to  7.2  km/sec. 
Since  the  velocities  are  substantially  different  and  since  there  is  no 
desire  to  imply  an  intercontinental  correlation,  the  present  study  will 
refer  to  this  horizon  as  the  Riel1  discontinuity  as  suggested  by  Professor 
D.  H.  Hall  of  Manitoba. 

The  reflection  labelled  M  was  of  erratic  quality  but  could  still 
be  correlated  over  parts  of  the  profile.  At  the  position  where  the  east- 
west  refraction  survey  was  recorded,  the  reflecting  horizon  corresponded 
in  depth  to  that  of  the  Mohorovi£id  discontinuity  as  determined  by  that 
survey.  Other  quasi-continuous  horizons  were  generally  of  poorer  quality 
and  some  of  these  exhibited  contrary  dip.  One  event  which  lies  0.9 
seconds  deeper  than  R  and  mirrors  the  latter's  structure  is  most  certainly 
of  multiple  origin.  This  is  shown  in  synthetic  seismogram  studies  (Clowes 
et  al.,  1968)  where  a  reverberation  between  the  free  surface  and  the  top 
of  the  Paleozoic  causes  a  first-order  multiple  to  occur  about  0.9  seconds 

^ouis  Riel,  a  colorful  western  Canadian  intellectual  and  pioneer 
of  French-Indian  ancestry  who  brought  the  province  of  Manitoba  into  union 
with  Canada,  became  a  member  of  parliament,  and  finally  led  a  revolt  which 
ended  his  career  on  a  scaffold. 
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after  the  primary  event . 

Intermingled  groups  of  many  closely  spaced  reflections  at  some 
positions  probably  show  the  effects  of  three-dimensional  structure  pro¬ 
jected  on  a  two-dimensional  cross  section.  Some  steeply  dipping  events 
could  be  due  to  energy  scattered  from  steep  structural  features  or  faults. 
Some  indication  of  three-dimensional  structure  was  obtained  from  a  2.4 
kilometer  east-west  profile  located  at  the  17-kilometer  distance  on  the 
diagram.  Good  continuous  reflections  which  corresponded  to  events  at 
similar  times  on  the  major  profile  were  recorded.  All  reflecting  hori¬ 
zons  showed  considerable  west  dip,  this  amounting  to  about  12  degrees 
for  the  Riel  discontinuity.  Such  results  emphasize  the  need  for  detailed 
investigations  and  imply  that  one  should  interpret  two-dimensional  sec¬ 
tions  with  caution  when  the  true  structure  probably  has  major  three- 
dimensional  complications. 

Structural  relief,  as  indicated  by  the  cross  section  of  Figure 
4.1,  is  of  considerable  magnitude.  The  horizon  interpreted  as  the  M 
discontinuity  varies  from  a  depth  of  47  to  38  kilometers.  At  the  same 
time,  the  Riel  discontinuity  changes  from  a  maximum  depth  of  37  kilometers 
to  a  minimum  of  26  kilometers.  The  latter  variation  represents  relief 
of  11  kilometers  over  a  horizontal  distance  of  only  50  kilometers.  Con¬ 
sidering  that  for  most  refraction  interpretations  in  stable  areas,  the 
Moho  is  assumed  to  be  planar  with  either  little  or  no  dip,  it  is  surpris¬ 
ing  and  enlightening  that  the  deeper  horizon  shows  comparable  relief  to 
the  shallower  one.  The  profile  suggests  that  the  isopach,  R  to  M,  is 
increasing  toward  the  north,  as  it  should  to  agree  with  depths  calculated 
by  Richards  and  Walker  (1959)  for  a  reversed  refraction  profile  centered 
about  60  kilometers  to  the  northwest.  Evidence  for  the  faults  could  be 
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seen  on  some  seismograms  where  good  quality  reflections  within  about  one 
second  of  each  other  and  having  similar  character  but  opposing  stepouts 
could  be  picked.  When  these  events  were  migrated  to  their  true  reflecting 
positions,  the  hypothesis  of  faulting  was  further  strengthened.  Although 
this  initial  interpretation  was  expected  to  be  somewhat  preliminary,  it 
has  stood  up  well  in  light  of  much  further  evidence. 

4.2  The  digitally  processed  seismograms 

The  dips  of  the  structures  indicated  on  the  two-dimensional  pro¬ 
file  presented  in  the  last  section  are  manifested  on  the  recorded  seismo¬ 
grams  as  reflection  events  with  varying  apparent  velocities.  Consequently 
the  form  of  digital  filter  described  in  Section  2.4 — the  velocity  filter 
which  selectively  enhances  reflections  within  a  specified  pass  range  of 
apparent  velocities— is  admirably  suited  for  delineating  any  possible 
coherent  reflected  energy  which  may  be  arriving  from  dipping  structures. 

Filters  of  this  type  have  been  applied  to  the  entire  suite  of 
digitized  reflection  seismograms  except  for  cases  where  one  of  the  two 
individual  recordings  necessary  for  the  filtering  process  was  not  avail¬ 
able.  In  genera^  the  beam-forming  filter  was  applied  with  three  differ¬ 
ent  pass  bands:  -11.7  to  -35.2  km/sec,  -35.2  to  +35.2  km/sec  and  +35.2 
to  +11.7  km/sec.  If  these  values  are  converted  to  time  stepout  per 
trace,  the  total  range  of  such  moveout  is  from  -25  to  +25  milliseconds 
per  trace  (subsurface  dips  of  about  -17°  to  0°  to  +17°  at  the  level  of 
the  Riel).  Whenever  plots  of  the  original  data  or  filtered  output  sug¬ 
gested  that  apparent  velocities  other  than  those  generally  used  might 
be  significant,  a  filter  with  different  characteristics  would  be  applied. 
In  a  few  cases  a  search  for  coherent  energy  with  very  low  apparent 
velocities  was  made,  but  this  provided  few  correlatable  events. 
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By  carefully  analysing  the  velocity  filtered  seismograms,  a 
complete  record  section  along  the  Lomond  and  Blackfoot  profiles  has  been 
obtainedo  Frequently  it  was  found  that  one  of  the  three  filtered  out¬ 
puts  showed  coherent  reflection  phases  with  greater  amplitudes  than 
either  of  the  remaining  two.  In  such  cases,  it  was  easy  to  choose  an 
appropriate  filtered  record  for  using  in  the  section.  However  quite 
often  it  was  found  that  one,  two  or  all  three  filtered  outputs  would 
have  some  reflections  with  good  signal-to-noise  ratios.  Usually  these 
would  be  at  different  traveltimes  and  it  was  necessary  to  piece  together 
a  composite  record  to  represent  all  the  reflection  events.  Of  course 
the  choice  of  which  filtered  seismogram  provided  the  best  data  was  not 
always  simple.  Complications  arising  from  the  simultaneous  arrival  of 
energy  representing  reflections  from  horizons  of  opposing  dips  could  be 
found  on  some  recordings.  To  further  confuse  the  interpreter^  the 
effects  of  three-dimensional  structure  producing  some  coherent  reflected 
energy  must  be  considered. 

In  Figure  4.2,  the  total  reflection  record  section  compiled  from 
the  velocity  filtered  seismograms  is  presented.  Because  of  the  large 
quantity  of  data,  the  profile  has  been  divided  into  six  sections  and  the 
illustration  continues  for  six  pages.  The  record  section  progresses 
from  south  to  north,  from  the  Lomond  to  the  Blackfoot  profile.  In  order 
to  provide  continuity,  the  last  record  on  each  page  of  the  section  is 
duplicated  as  the  first  record  of  the  succeeding  page  (except  for  4.2B 
to  4.2C  where  there  was  one  set  of  seismograms  unavailable).  As  an  aid 
for  comparison  with  the  cross  section  of  Figure  4.1,  the  caption  for 
the  present  diagram  lists  the  approximate  distance  ranges  according  to 
the  scale  of  Figure  4.1  for  each  part  of  the  illustration.  Blank  spaces 
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Figure  4.2  Velocity  filtered  reflection  record  section  along 
a  90 -kilometer  -profile .  On  the  labels  to  the 
left  of  each  seismogram s  the  first  two  numbers 
refer  to  the  records  which  were  filtered.  The 
numbers  following  "VF"  give  the  range  of  apparent 
velocities  for  the  pass  band  of  the  velocity  filter . 

If  the  results  from  more  than  one  pass  band  have 
been  considered this  is  indicated  at  the  appropriate 
position  above  the  seismogram „ 

The  illustration  continues  for  six  pages  and  the 
appropriate  part  of  the  cross-section  of  Figure  4°  l  to 
which  each  page  refers  is  shown  below  according  to 
the  distance  scale  of  Figure  4. 1. 


4.2A  10  to  24  kilometers . 

4.2B  22  to  35  kilometers. 

4.2C  36  to  50  kilometers . 

4. 2D  50  to  71  kilometers . 

4„2E  69  to  87  kilometers . 

4.2F  85  to  102  kilometers . 


Note  that  Figure  4.1  is  a  migrated  structure  section 
while  the  foregoing  distances  refer  to  the  surface 
recording  locations . 


' 


.i  ‘ 


i  •!  tc  •'  f.  t/tf  vr  »  t  ^  .**  ^  ' 


RECORD  09  1966  7,8  VF  -35  TO  35  K/S  5.6  VF  -35  TO  35  K/S  3.4  VF  -35  TO  35  K/S 


8.00  e.90 


19.90  19. so  is.  M  is.® 


1 

>\/jV-Ay\AA /'^/W^/vv\ll/i/l/^^  >*  \x 

/VvV^vrVl^^  7\/VVw\JWa/vv  “v\mt  \* 

»A.«.  .  l.i  K„.  .11.1  «  .  .A  »  »  fl  A..  .  .  A  *..»,  .  .  »  ».  J|yy\/\/AAA/K/UA/VA^'  \V\/''V'/V< 


VF  -9,-18  KM/S 


AVvWV^wVAyVyVy^YxA/^ 
*'^/\A/\^\^~AA/'vA^YAA^V^^ 
X'^V\MAat~^a~/W'1/'yAAWI/\/\cY^'WV^^ 


/yYV'M/vVWV\^/vvMMAAM/VVvv"v\^^ 
Jy\l\J\j\W^wvVVJ^\/^[\/foJ\J\/J\JK/^\/\hKt\/^^^^ 

— v^A/\/T_^AA^^Vv^VVV^< 


8.00  8.90  8.80  9.20  9.60  10.00  10.90  10.80  11.20  11.60  12.00  12.90  12.80  13.20  13.80  19.00  19.90  19.80  IS.20  IS. 80 

TIME  (SECONDS) 


29.30  VF  -35.35  KM/S  27,28  VF  -12.-35  K/S  25.26  VF  -12.-35  K/S  23.24  VF  -12.-35  K/S  21.22  VF  -12,-35  K/S  19,20  VF  -12,-35  K/S  17.18  V F  -12.-35  K/S  15.16  VF  35  -  12  K/S 


8.00  8.40 


9.20  9.80  10.00  10.110  10.80  11.20 


TIME  (SECONDS) 


12.40  12.80  13.20  13.80 


4.80  IS. 20  IS. 80 


8 


10 


14 


A 


M 


^vT\|V'AjW  — x 

am Hwim 


L  \M/w 


i 

irvMn 


VW  vWi/wMN 


^Yv\A\P^Wm/Wvvv«« 
\/MA\^/vvWV/\rtM-x 

..  - « ll4lVvWWrv^^AA/  V\/\^A aA^vAAAvA^< 


'^WVMAA'"  ^ 


1 

i\L 


VM/Wm/I/ 


uaLaVaTmaA 

vviHvviur  m 


\|  y\/Vv\/W'V\^n/WV^\/^/VWA^x 


\m/W\A^]Wy^^ 


v  vy  v 


Mf 


/vvyvvu 


A  A  A|\/\/\^W\/WVWvv< 


^  ^  ^  (nVWW\AAAA|VAAA/W\/^v^vv^ 


^AA/VW^/vywv/WWVV^^ 


xyA^yvWv\ft^vN/VlAAA/\A/'. 


\AA^v^^A<v\/^AAA7\AA/\/yVVVAA^\^ANvW\Av^^-~^< 
AvvAA^AA/VVArtAA^Af^/^MAF^v^^** 
V\MAAV \aTHAN]^I\J{a^^ 


8  10  12  14  16 

I _ . _ . _ - _ . _ I _ . _ _ _ l-  . _ _ _ . _ I _ . _ . _ . _ . _ I 

8.00  9.40  8.80  9.20  9.80  10,00  10.40  10.80  11.20  11.80  12.00  12.40  12.80  13.20  13.80  14.00  14.40  14.80  >5.20  IS. 80 

TIME  (SECONDS) 


,-y 


mmam mt* ' 
mSm.^ 


47.143  VF  -12.-35  K/S  45,46  VF  -12,-35  K/S  43,44  VF  -12.-35  K/S  41,42  VF  -12.-35  K/S  39,40  VF  -12,-35  K/S  37.38  VF  -35,35  KM/S  35.36  VF  -35,35  KM/S  33,34  VF  -12,-35  K/S 


8.40  8.80  9.21 


4.40  14.80  IS. 20  15.80 


8 


10 


12 


14 


16 


|  R  l  M 

lil  l\f lA  M 1  AIlllIW^N/w'A/y\rA/^Wl/v\/wMAAA/'/A^^^^ 


''\/\/VV\A'-/Wv'''VVv-^^ 

WVVv'AAAAA/\AA/MA^wvwwWv/wwvyV\A-»< 

V\AAAA/\MMMA/Wv - aAa-vuiAm< 


14.80  IS. 20  IS. 80 


' 


,  ,  „  JJpw  ■  Mmiwfffc 

timrvim^,Smw 


Hwf fwfaKmftm  ’*  pjv&  ^jwfW 
MM0 ' m  '*if{ fM&fj ( 


mmMM 


36.35  VF  -12, -35KM/S  H2.H1  VF  -35.35  KM/S  HH.H3  VF  -12.-3SKM/S  H6.H5  VF  -12.-3SKM/5  S0.H9  VF  -12.-35KM/S  S2.51  VF  -35.35  KH/S  56.55  VF  -11.-3SKM/S 


J  ■  J  . .  ,  ,  .  .  I 


I  Fit  (SCXONOS) 


IS. 2D  IS. ID 


14 


i.OO  8. no 


TIME  ISECONOS) 


OO  I2.no  12.80  13.20  13.80 


AnM^ 

xv^y  \ !\[ v^X^AA/y^ATy^v/  V\ 


vw'a~A^M/x 


\Z\>\^yy\i^Aj^v,'s^s*AS — virfi 

</yUv/'^Y'W'''H/'A/Ww^WVV^ 


mmim « 


®s*j 

1  MM  $*»{ 

#w ! 


»*  Si ,k«  y4H* 

V* 


02.01  VF  -12. -3SKM/S  09.03  VF  -12.-35KM/5  10,09  VF  -35.35  KM/S  12. 11  VF  -35.35  KM/S  19. 13  VF  -35,35  KM/5  16,15  VF  35.12  KM/S  18. 17  VF  35. 12  KM/S  20. 19  VF  -35.35  KM/S 


Aw  *  * » 

. ,.  Ji .  mmm WSkm**  f* 

f^Wn  w*f*M*Y  #Prt;A  *  '•  rWSf 


|»P  I 

P*->  T.?»  ‘IS 

■??.#  'Mf  '«. 


95 


on  the  seismograms  represent  time  ranges  for  which  there  were  instrumental 
difficulties  or  large  amplitude  low  frequency  waves.  The  few  cases  where 
the  latter  interference  effects  have  been  left  in  the  illustration  are 
obvious  from  their  low  frequency  and  unusual  character.  In  one  instance 
(Record  09,  1966  of  Figure  4.2A)  where  only  one  of  the  two  seismograms 
necessary  for  filtering  was  available,  the  original  unfiltered  data  has 
been  substituted  for  the  sake  of  continuity. 

While  making  a  correlation  across  the  fan-pass  filtered  seismo¬ 
grams,  it  must  be  recalled  that  the  filtering  process  reduces  the  number 
of  traces  from  twelve  to  five.  If  the  original  traces  were  numbered  1  to 
12  to  reference  their  surface  position,  the  filtered  channels  would  be 
labelled  4.5  to  8.5  in  which  the  trace  at  4.5  is  a  weighted  composite  of 
traces  1  to  8,  the  one  at  5.5  of  traces  2  to  9,  etc.  On  the  basis  of 
the  horizontal  separation  between  station  locations,  the  spacing  between 
adjacent  seismograms  should  be  2.4  times  greater  than  that  depicted  to 
present  a  true  scaled  section.  This  fact  should  be  noted  when  correla¬ 
tions  are  being  made.  In  addition, consideration  should  also  be  given  to 
the  results  concerning  transition  zones  which  were  discussed  in  Chapter  3. 
It  is  possible  that  the  largest  amplitude  peak  in  a  given  reflected  sig¬ 
nal  sequence  could  shift  from  record  to  record  depending  on  the  thickness 
of  the  transition  zone.  However, since  this  must  occur  within  a  depth 
interval  of  less  than  one  kilometer,  for  any  reasonable  velocity  values 
chosen  the  oscillations  of  the  pulse  series  should  be  negligible  within 
less  than  0.4  seconds. 

The  reflection,  R,  from  the  Riel  discontinuity  is  usually  dis¬ 
tinguished  as  the  most  prominent  event  on  the  seismograms  and  as  such, 
it  could  be  correlated  most  easily  across  the  entire  section.  Two  peaks 
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of  the  pulse  sequence  associated  with  R  have  been  darkened  to  more  clearly 
show  the  correlation.  In  many  cases,  these  are  followed  by  another  large 
amplitude  pulse  about  one  second  later  and  this  is  most  likely  of  multiple 
origin.  With  only  a  few  exceptions,  a  reflection,  M,  interpreted  as  being 
from  the  Mohorovi<5id  discontinuity  has  also  been  correlated  across  the 
suite  of  records.  While  in  many  cases  signal-to-noise  ratios  were  only 
slightly  greater  than  one,  the  correlation  is  aided  by  considering  other 
sequences  of  pulses  and  their  variation  from  record  to  record.  Peaks 
from  this  event  have  also  been  darkened  for  emphasis.  For  both  R  and  M, 
correlations  at  the  northern  end  of  the  Blackfoot  profile  (Figures  4.2E 
and  4.2F)  were  more  ambiguous  and  a  definitive  interpretation  is  difficult. 

One  of  the  most  obvious  results  to  be  derived  from  the  record 
section  is  the  large  number  of  events  which  have  been  revealed  by  the 
velocity  filtering  process.  Some  of  these  could  be  correlated  over  much 
of  the  profile.  Two  of  the  more  continuous  sequences,  one  in  the  inter¬ 
val  between  R  and  M  and  one  about  1.5  seconds  prior  to  R,  have  been 
designated  by  single-headed  arrows.  In  most  cases  the  existence  of  these 
reflections  is  quite  clear  and  in  a  few  instances  an  excellent  quality 
reflected  pulse  has  been  delineated.  On  some  seismograms  for  which  the 
dynamic  range  of  the  instrumentation  was  not  exceeded  and  surface  waves 
were  negligible,  much  shallower  reflections  were  detected.  In  two  such 
cases  (Records  32,31  and  28,27  of  Figure  4.2E),  a  good  quality  reflec¬ 
tion  was  obtained  just  five  seconds  after  the  shot  instant.  On  the  part 
of  the  record  section  where  the  structure  reaches  its  maximum  height 
(Figures  4. 2D  and  E) ,  indication  of  a  continuously  correlatable  pulse 
can  be  detected  at  about  13.8  seconds,  nearly  two  seconds  later  than  the 
reflection  attributed  to  the  M  discontinuity.  Over  the  ten  seismograms 
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on  which  it  is  quite  clear,  the  pulse  has  been  marked  with  a  double-headed 
arrow.  Further  evidence  for  reflected  energy  from  below  the  Moho  is 
given  on  Records  05,06  and  07,08  (Figure  A. 2k)  near  the  southern  end  of 
the  profile.  A  good  amplitude  event  is  seen  at  a  two-way  reflection 
time  of  about  16.8  seconds.  Beloussov  et  al .  (1962)  and  Khalevin  et  al. 
(1966)  have  reported  similar  observations  from  the  deep  seismic  sounding 
program  in  the  Soviet  Union  where  good  reflections  as  late  as  18  seconds 
have  been  observed.  It  is  possible  that  the  late-arriving,  steeply 
dipping  event  could  be  associated  with  diffracted  energy  from  the  fault 
postulated  at  the  southern  end  of  the  Lomond  profile.  On  the  other  hand, 
some  of  the  relfections  from  low  dipping  interfaces  may  be  from  skills 
within  the  upper  mantle. 

The  velocity  filtered  data  have  confirmed  the  gross  structural 
features  of  the  cross  section  derived  in  Section  4.1  and  have  provided 
additional  details.  For  example,  the  R  and  M  horizons  have  been  more 
reliably  correlated  across  the  profile  and  indicate  that  the  structure 
of  the  horizons  is  less  smooth  than  depicted.  The  correlation  of  other 
reflections  across  much  of  the  profile  provides  evidence  of  more  extended 
layering  than  is  usually  considered  within  the  deep  crust.  These  horizons 
also  exhibit  similar  structural  variations.  Some  support  has  been  ob¬ 
tained  for  the  possible  existence  of  reflecting  interfaces  below  the  M 
discontinuity.  As  well,  some  evidence  is  given  for  block-faulting  in 
contrast  to  the  single  fault  plane  interpreted  in  the  original  cross 
section.  For  clarity.  Figure  4.3  illustrates  three  filtered  seismograms 
which  show  the  good  correlation  across  the  central  fault  of  the  section. 
The  offset  represents  a  time  difference  of  about  one  second.  This  inter¬ 
pretation  is  consistent  with  refraction  data  which  is  discussed  in  the 
following  section. 
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Figure  4. 
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Even  without  band-pass  filtering,  some  of  the  fan-pass  filtered 
seismograms  still  have  an  oscillatory  appearance.  For  such  cases  the 
reflection  events  could  possibly  be  more  clearly  distinguished  by 
utilization  of  de reverberation  filters.  The  application  of  this  sug¬ 
gestion  would  be  a  difficult  procedure  as  a  fairly  complete  knowledge 
of  the  reflection-producing  capabilities  of  the  sedimentary  layering 
would  be  required.  However  the  present  section  has  shown  that  signi¬ 
ficant  improvement  in  record  quality  is  rendered  by  the  application  of 
just  the  velocity  filter  to  the  observed  data. 

4.3  The  refraction  seismograms 

Interpretation  of  the  reflection  observations  has  shown  that 
large-sized  structures  exist  on  the  Mohorovicid  discontinuity,  within 
the  lower  crust  and  probably  extend  into  the  upper  crust.  This  implies 
that  it  should  be  possible  to  detect  such  features  with  a  detailed 
refraction  study.  For  the  broadside  arc.  Profile  C,  of  the  refraction 
survey,  the  positions  of  the  seismic  detectors  were  selected  to  most 
clearly  detail  those  areas  where  faulting  was  suspected  (see  Figure  1.1 
for  recording  station  locations  and  profile  designations) .  On  the 
twelve  useable  seismograms,  two  prominent  refraction  events  could  be 
clearly  seen.  The  earlier  of  these  arrivals  corresponded  to  Pn,  a 
compressional  wave  travelling  in  the  upper  mantle.  The  second  arrival, 
designated  P  ,  was  interpreted  as  a  head  wave  from  within  the  Precambrian 

o 

basement  complex  on  the  basis  of  a  reversed  refraction  profile  previously 
obtained  in  southern  Alberta. 

The  in-line  profiles,  A  and  B,  were  chosen  to  observe  differences 
in  phase  velocity  and  intercept  times  over  the  deepest  part  of  the 
structure  and  on  the  highest  part  to  the  north.  Figure  4.4  shows  the 
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Figure  4.4  Traveltime  curves  for  Pg  and  Pn  from  profiles 
A  and  B  (Figure  l.l) .  For  Profile  A  recorded 
along  the  center  of  the  deep  part  of  the  structure 3 
the  arrivals  are  about  one  second  later  than  for 
Profile  B  recorded  along  the  structural  high. 
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travel time  curves  for  Pg  and  Pn.  It  is  noted  that  there  is  a  distinct 
difference  in  the  intercept  time  of  the  two  profiles,  although  the  phase 
velocities  are  only  slightly  different.  The  high  apparent  velocity 
(8.67  km/sec)  for  Pn  is  probably  indicative  of  dip  since  the  reversed 
profiles  have  been  interpreted  as  having  a  velocity  of  8.3  km/sec.  In 
Figure  4.5,  the  residual  Pg  and  Pn  traveltimes  for  Profile  C  are  plotted. 
Nearly  a  one-second  difference  in  Pn  arrival  times  and  about  1.3  seconds 
difference  in  Pg  arrival  times  occurs  over  a  north-south  horizontal 
distance  of  less  than  ten  kilometers.  Stations  C3  and  A1  are  recording 
arrivals  from  the  upper  and  lower  sides,  respectively,  of  the  postulated 
fault  in  the  center  of  Figure  4.1.  Although  the  amount  of  structural 
variation  to  be  interpreted  from  such  a  change  is  critically  dependent 
on  the  crustal  model  chosen,  the  general  result  certainly  confirms  the 
existence  of  a  steeply  dipping  major  fault  in  the  crust.  The  similarity 
in  reduced  traveltime  curves  for  Pg  and  Pn  indicates  that  the  fault  does 
indeed  extend  throughout  the  entire  crustal  section  of  Precambrian  age. 

4.4  The  gravity  and  magnetic  observations 

The  gravity  survey  carried  out  by  the  Geophysics  Division  recorded 
values  from  a  Worden  Model  XP  0  gravimeter.  Normal  field  procedures  were 
followed  and  the  data  were  reduced  to  Bouguer  gravity  anomalies  following 
the  method  adopted  by  the  Gravity  Division  of  the  Dominion  Observatory. 
This  survey  was  then  combined  with  the  Bouguer  gravity  data  provided  by 
the  Gravity  Division.  Their  compilation  is  comprised  of  the  base  stations 
established  by  the  Division  plus  values  submitted  by  Gulf  Canada,  Limited 
and  Chevron  Standard,  Limited,  two  oil  companies  that  have  worked  in  the 
area.  In  cases  where  university  and  government  stations  coincided,  it 
was  found  that  Bouguer  values  from  the  former  agreed  within  a  fraction 
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Figure  4.5  Residual  Pg  and  Pn  traveltime  curves  for 
Profile  C  (Figure  l.l) .  The  horizontal 
distance  between  stations  C3  and  A1  is  less 


than  ten  kilometers. 
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of  one  milligal  with  those  listed  by  the  government.  All  gravity  sta¬ 
tion  locations  which  were  used  in  compiling  the  map  are  indicated  on 
Figure  4.6,  with  approximately  one- third  of  these  occupied  by  the 
university. 

The  complete  Bouguer  gravity  anomaly  map  is  presented  in  Figure 
4.7.  To  the  northeast  of  Brooks,  there  is  a  small  but  prominent  posi¬ 
tive  anomaly,  called  the  Princess  high,  which  is  associated  with  a 
local  structural  high  on  the  Precambrian  surface.  This  feature  has  been 
penetrated  by  a  well  and,  from  a  small  sample  of  the  material,  it  is 
known  that  the  density  is  high — about  3.0  gm/cm3  (R.M.  Burwash,  personal 
communication).  However,  the  seismic  reflection  profiles  are  a  consider 
able  distance  west  of  this  anomaly  and  are  probably  unaffected  by  it. 

Two  features  relative  to  these  profiles  are  important  to  note.  First, 
in  the  region  of  the  Lomond  profile,  there  is  a  pronounced  east-west 
trending  Bouguer  gravity  low  which  cuts  across  most  of  southern  Alberta, 
although  it  is  distorted  by  the  local  high  discussed  above.  To  the 
west,  the  course  of  the  anomaly  is  intercepted  by  the  gravity  effects 
of  the  Rocky  Mountains,  which  trend  nearly  north-south.  The  second 
major  feature  of  interest  for  this  project  is  the  broad  gravity  high 
centered  on  the  northern  end  of  the  Blackfoot  profile.  Its  large 
expanse  and  considerable  magnitude,  suggestive  of  deep-seated  structural 
variations,  helped  determine  the  location  of  the  latter  profile. 

For  the  ground  magnetometer  survey  run  by  the  Geophysics  Division 
values  of  the  total  magnetic  field  intensity  were  read  from  a  Barringer 
Model  GM-102A  nuclear  precession  magnetometer  with  an  accuracy  of  about 
ten  gammas.  The  data  recorded  in  the  field  were  adjusted  to  minimize 
any  differences  in  readings  at  repeated  stations.  By  computing 
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Figure  4.6  Gravity  station  locations  for  the  map  of 

Figure  4.7.  The  heavy  lines  show  the  loca¬ 
tion  of  the  Lomond  and  Black  foot  -profiles 
and  the  east-west  expanding  spread  profile. 
The  dashed  lines  outline  the  assumed  bounda¬ 
ries  of  the  ancient  rift  valley  (Section  4.5) 
as  traced  from  magnetic  and  gravity  trends. 
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Figure  4.7  Bouguer  gravity  anomaly  map  for  southern 
Alberta  and  southeastern  British  Columbia. 
The  contour  interval  is  5  milligals . 


106 


the  magnetic  field  due  to  dipole  and  higher  order  effects  by  means  of 
a  seventh-degree  spherical  harmonic  analysis,  then  subtracting  these 
values  from  the  observed  ones,  a  residual  total  magnetic  intensity  map 
was  derived  (Figure  4.8).  The  small  dots  represent  the  locations  at 
which  magnetometer  readings  were  made.  This  residual  magnetic  map  is 
also  characterized  by  strong  east-west  trending  lineaments  across  the 
region  of  the  two  profiles.  Although  it  is  centered  slightly  further 
south,  the  large  negative  anomaly  corresponds  with  a  similar  feature  on 
the  Bouguer  gravity  map.  Of  particular  interest  is  the  continuation  of 
the  east-west  trends  underneath  the  Rocky  Mountains  which  lie  almost 
perpendicular  to  these  trends.  It  shows  that  the  magnetic  effects  of 
the  more  recent  mountain  system  are  negligible  and  causes  of  the  resid¬ 
ual  anomalies  must  be  associated  with  deep-seated  phenomena. 

On  the  structure  cross  section  of  Figure  4.1,  the  deepest  part 
of  the  features  corresponds  with  both  a  gravity  and  magnetic  low,  while 
the  horst  to  the  north  coincides  with  both  a  gravity  and  magnetic  high 
(upper  part  of  Figure  4.1).  Figure  4.9  shows  the  results  of  gravity 
calculations  based  on  a  two-dimensional  generalization  of  the  seismic 
section.  The  observed  Bouguer  gravity  (solid  curve)  and  residual 
magnetic  (dashed  curve)  profiles  are  also  given.  Velocities  in  the 
post-Precambrian  sediments  were  obtained  from  well  log  data  while  those 
in  the  deeper  crust  were  taken  from  the  results  of  the  reversed  refrac¬ 
tion  surveys  mentioned  previously.  Corresponding  rock  densities  were 
selected  from  these  velocities  on  the  basis  of  a  model  proposed  by 
Nafe  and  Drake  (1959).  The  small  circles  on  the  observed  Bouguer 
profile  show  the  fit  of  the  proposed  crustal  section.  By  seismic  means, 
the  existence  of  the  fault  in  the  central  part  of  the  section  has  been 
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Figure  4.8  Residual  total  magnetic  field  intensity  map 

for  southern  Alberta  and  southeastern  British 


Columbia . 
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Figure  4.9  Generalized  crustal  section  across  southern 
Alberta .  Dashed  curve:  residual  total 
magnetic  field  intensity .  Solid  curve:  Bouguer 
gravity  anomaly .  The  small  circles  on  the 
latter  shew  the  fit  of  the  proposed  crustal 
model  depicted  in  the  lower  half  of  the  figure. 
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confirmed,  but  the  postulated  fault  on  the  southern  end  of  the  profile 
was  based  on  ambiguous  reflection  data.  The  gravity  calculation  shows 
that  this  fault  is  required  to  produce  the  steep  gradient  observed 
directly  above  it.  The  gravity  and  seismic  data  are  best  fitted  if  the 
two  fractures  are  taken  to  be  steeply  dipping  faults  that  extend  to  the 
top  of  the  Precambrian  section  at  a  depth  of  about  2.5  kilometers.  On 
the  basis  of  several  wells  drilled  to  the  basement  complex,  it  is  known 
that  no  major  faults  continue  into  the  Paleozoic  sediments.  From  these 
considerations  it  seems  reasonable  to  describe  the  deepest  part  of  the 
structure  as  a  graben  within  the  ancient  crust.  This  graben  must  be 
filled  with  low  density  Precambrian  rocks  to  satisfy  the  gravity  low  and 
an  upper  crustal  refracting  horizon  having  a  low  apparent  velocity.  At 
the  level  of  the  Mohorovicid  discontinuity,  mass  variations  do  not  have 
a  great  effect  on  gravity  models  so  the  premise  of  low  density  fill  is 
necessary  to  account  for  the  pronounced  gravity  low. 

Garland  and  Burwash  (1959)  found  a  correlation  between  the 
lithology  of  the  Precambrian  basement  and  gravity  observations  in  central 
Alberta.  From  the  data  just  presented,  it  appears  that  a  similar  cor¬ 
relation  exists  in  southern  Alberta  between  lithology  as  suggested  by 
seismic  velocities  and  both  gravity  and  magnetic  data.  The  seismic  sec¬ 
tion  indicates  that  such  correlations  are  directly  due  to  deep  crustal 
structure  in  which  different  blocks  have  been  raised  or  lowered  by  as 
much  as  ten  kilometers.  In  the  present  case,  there  appears  to  be  a 
graben  of  Precambrian  age  which  was  filled  with  lighter  density  rocks 
having  a  lower  proportion  of  magnetically  susceptible  mineralization. 

Thus  there  exists  a  direct  relationship  between  the  gravity  and  magnetic 
profiles . 
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4o5  A  buried  Precambrian  rift 

The  most  prominent  feature  of  the  generalized  crustal  profile 
presented  in  the  last  section  is  a  graben  coinciding  with  gravity  and 
magnetic  lows  which  form  part  of  an  east-west  trending  lineation.  From 
this  correlation,  the  graben  appears  to  be  a  major  east-west  structure 
extending  at  least  450  kilometers  across  southern  Alberta  and  passing 
at  right  angles  to  the  strike  of  the  more  recent  Rocky  Mountain  system. 
The  latter  are  underlaid  by  a  peneplaned  and  gently  westward  sloping 
crystalline  basement  (Bally  et  al.,  1966)  which  has  pot as si urn- argon 
dates  of  1.5  to  2.0  billion  years  (Burwash  et  al.,  1962)  and  forms  a 
portion  of  the  Churchill  geological  province  (Kanasewich,  1965).  From 
the  evidence  of  its  great  elongation  on  the  gravity  and  magnetic  maps, 
the  graben  may  best  be  described  as  a  rift.  Since  post-Middle  Cambrian 
sediments  are  not  known  to  be  affected  by  the  faulting,  it  seems  reason¬ 
able  to  conclude  that  this  structure  represents  a  buried  Precambrian 
rift.  The  probable  location  of  this  feature  as  traced  with  the  aid  of 
the  gravity  and  magnetic  lineaments  is  sketched  in  Figure  4.6. 

While  the  existence  of  Precambrian  rifting  has  not  previously 
been  established,  there  are  good  examples  of  more  recently  developed 
continental  fractures  such  as  the  East  African,  Rhinegraben  and  Lake 
Baikal  rift  systems.  The  rift  valleys  of  East  Africa  have  been  exten¬ 
sively  discussed  (Gregory,  1896  and  1921),  but  only  limited  geophysical 
studies  have  been  made  (see  Irvine,  1966).  Both  the  Rhinegraben  (see 
Rothe  and  Sauer,  1967)  and  the  Lake  Baikal  (see  Bulmasov,  1960;  Zorin, 
1966;  Reisner,  1966;  and  Ladynin,  1966)  rift  systems  have  undergone  more 
extensive  geophysical  analyses.  Thus  some  comparisons  of  the  recent 
rifts  can  be  made  with  the  Precambrian  rift  in  western  Canada. 
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Freund  (1966)  notes  that  most  continental  and  oceanic  rifts  have  widths 
ranging  from  30  to  70  kilometers,  although  there  are  exceptions.  The 
width  of  the  rift  whose  outlines  are  traced  in  Figure  4.6  averages  about 
40  kilometers.  At  Lake  Baikal,  the  maximum  overall  depth  of  the  rift, 
from  the  highest  mountain  crest  to  the  surface  of  the  sediment  fill,  is 
greater  than  4.1  kilometers  (Hope,  1967).  A  similar  value,  4.4.  kilo¬ 
meters,  is  established  by  lilies  (1967)  for  the  total  throw  of  the 
faulting  in  the  Rhinegraben.  In  the  structure  sections  of  Figures  4.1 
and  4.9,  the  offsets  along  the  fault  amount  to  about  3  to  5  kilometers. 
With  regard  to  the  sedimentary  fill,  lilies  (loc.  cit.)  estimates  a 
maximum  depth  extent  of  3.4  kilometers,  while  Florensov  (1966)  places 
the  total  thickness  of  sediments  in  the  Baikal  rift  zone  at  about  5  to 
6  kilometers.  On  the  basis  of  gravity  and  refraction  surveys,  the 
graben  fill  of  Figure  4.9  varies  in  thickness  from  about  3  to  7  kilo¬ 
meters.  From  an  analysis  of  deep  reflection  seismograms,  Demnati  and 
Dohr  (1965)  concluded  that  the  Conrad  discontinuity  was  depressed 
beneath  the  Upper  Rhine  valley  and  suggested  about  5  kilometers  of 
relief.  From  geophysical  work  on  the  delta  of  the  Selenga  River  and 
surface  geological  mapping,  Florensov  (loc.  cit.)  places  the  total 
amplitude  of  crustal  displacement  at  about  5  to  7  kilometers.  These 
estimates  are  not  dissimilar  from  the  10  kilometers  of  total  displace¬ 
ment  obtained  for  the  rift  system  in  southern  Alberta  on  much  stronger 
geophysical  evidence.  Thus  it  is  seen  that  in  terms  of  physical  magni¬ 
tudes,  the  characteristics  of  the  Precambrian  rift  bear  a  remarkable 
similarity  to  those  determined  for  the  Rhinegraben  and  Lake  Baikal  rift 
zones.  It  is  unfortunate  that  more  comprehensive  geophysical  analyses, 
including  detailed  seismic  investigations,  have  not  been  made  in  the 
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region  of  the  East  African,  Rhinegraben  and  Lake  Baikal  rifts  in  order 
that  a  better  comparison  of  such  physical  characteristics  could  be  made. 

However,  both  gravity  and  magnetic  features  of  the  recent  rifts 
have  a  similar  expression  over  the  ancient  one.  In  all  three  cases,  a 
pronounced  gravity  low  parallels  the  central  region  of  the  rifting. 

Closs  and  Plaumann  (1967)  give  a  relative  magnitude  of  -20  to  -50  milli- 
gals  for  the  depression  in  the  Rhinegraben,  values  which  fit  well  with 
the  data  given  in  Figure  4.7.  From  a  consideration  of  two-dimensional 
models,  they  conclude  that  the  Bouguer  anomaly  across  the  graben  can  be 
fully  explained  on  the  basis  of  the  sedimentary  fill  in  the  upper  crust, 
and  it  is  not  necessary  to  invoke  causes  of  the  anomaly  being  situated 
in  the  lower  crust  or  upper  mantle.  Zorin  (1966)  reaches  a  similar 
conclusion  for  the  observed  anomalies  at  Lake  Baikal  and  additionally 
shows  that  the  edges  of  the  rift  are  bounded  by  strong  positive  anomalies. 
All  these  characteristics  were  noted  from  the  gravity  study  in  the  region 
of  the  Precambrian  rift.  Just  as  with  the  gravitational  field,  the 
magnetic  field  above  the  Rhinegraben  and  Lake  Baikal  is  characterized  by 
linear  trending  anomalies,  parallel  to  the  strike  of  the  rifts.  Bosum 
and  Hahn  (1967)  find  that  in  the  graben  negative  anomalies  exist  and  are 
bounded  by  greater  amplitude  positive  values  of  the  magnetic  field.  Over 
the  Baikal  folded  region,  Bulmasov  (1960)  finds  a  regional  field  which 
is  negative  although  strong  linear  positive  anomalies  stand  out  against 
this.  These,  he  argues,  are  caused  by  intrusions  of  basic  rocks  along 
faults  located  below  the  surface  (depth-faults).  A  strong  positive 
magnetic  field  exists  on  one  side  of  the  rift  zone.  On  the  residual 
magnetic  field  map  of  Figure  4.8,  negative  anomalies  are  evident  above 
the  Precambrian  rift  with  positive  residuals  bounding  it  on  either  side. 
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The  entire  map  is  dominated  by  the  east-west  trending  lineaments  associ¬ 
ated  with  the  rift.  Thus  it  is  found  that  the  prevalent  gravity  and 
magnetic  features  of  the  Rhinegraben  and  Lake  Baikal  bear  a  striking 
resemblance  to  those  associated  with  the  rift  in  southern  Alberta. 

That  a  rift  in  the  geologically  ancient  crustal  section  has  been 
discovered  possesses  some  significance,  for  DeSitter  (1959)  notes  that 
"the  great  rift  valleys  of  the  earth  are  all  relatively  young,  or 
rejuvenated  structures  ....  The  reason  for  the  fact  that  no  large 
rift  valleys  of  Hercynian  or  earlier  age  are  known  remains  obscure." 

The  recognition  of  these  recent  rift  valleys  is  generally  based  on  their 
spectacular  geomorphological  features  and  not  on  geologic  structure. 
However,  there  is  now  some  geological  evidence  for  the  existence  of 
Precambrian  rifting  in  East  Africa  (McConnell,  1967).  In  the  present 
research,  the  buried  Precambrian  rift  has  been  identified  on  the  basis 
of  geophysical  studies.  This  has  shown  evidence  for  large  amounts  of 
relief  with  steeply  dipping  surfaces  on  the  Precambrian  terrain.  If 
it  is  assumed  that  the  Precambrian  terrain  exposed  in  many  shield  areas 
has  undergone  similar  amounts  of  warping  during  the  long  period  of 
tectonic  evolution,  then  the  evidence  from  this  study  indicates  that  in 
some  places  on  an  ancient  peneplaned  craton,  it  is  possible  to  observe 
directly  the  rocks  ten  kilometers  deeper  on  one  section  than  on  an 
adj  acent  one . 
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CHAPTER  5 

CONCLUSION 

This  thesis  has  demonstrated  that,  by  careful  design  of  field 
procedures  including  the  use  of  patterns  of  holes  and  arrays  of  seis¬ 
mometers,  near-vertical  incidence  reflected  energy  from  within  and  at 
the  base  of  the  earth's  crust  can  be  successfully  recorded.  Velocity 
filtering  techniques  proved  to  be  an  effective  method  of  enhancing 
individual  reflection  events  and  thus  made  correlations  over  longer 
distances  more  reliable.  Not  only  do  the  observed  reflections  provide 
detailed  knowledge  concerning  structural  variations,  but  inherent  in 
the  reflected  events  is  information  relating  to  the  intrinsic  proper¬ 
ties  of  the  material  through  which  the  seismic  waves  have  propagated. 

The  interpretation  of  the  recorded  seismograms  has  revealed  the 
presence  of  large-scale  structural  features  and  steeply  dipping  faults 
within  the  lower  crustal  section  of  southern  Alberta.  By  combining 
gravity  and  magnetic  data  with  the  seismic  results,  a  major  east-west 
rift  of  Precambrian  origin  has  been  discovered.  Under  2.5  kilometers 
of  younger  sediments,  the  feature  extends  for  over  450  kilometers  and 
has  11  kilometers  of  relief  from  the  highest  to  the  lowest  part. 
Geomorpho logical  features  and  geophysical  characteristics  of  this 
ancient  rift  bear  a  remarkable  similarity  to  those  of  more  recently 
developed  continental  rift  valleys. 

On  the  basis  of  the  discovery  of  the  Precambrian  rift  described 
in  this  study,  and  from  other  data,  Kanasewich  (1968)  suggested  a 
possible  origin  for  the  lead- zinc  field  near  Kimberley,  British  Columbia 
and  the  Coeur  d'Alene  mining  district  of  Idaho.  He  postulated  that 
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these  ore  bodies  were  deposited  in  the  Precambrian  rift  under  conditions 
closely  related  to  those  prevailing  in  the  present  Red  Sea  rift,  where 
three  pools  of  hot  acidic  brines  with  extremely  high  concentrations  of 
heavy-metals  have  been  found. 

The  importance  of  crustal  effects  on  seismic  arrivals  recorded 
at  the  LASA  array  in  Montana  has  been  demonstrated  by  Mack  (1969) .  He 
considered  a  simplified  version  of  the  crustal  section  in  southern 
Alberta  (Clowes  et  al.,  1968)  to  show  that  such  deep  crustal  relief 
could  cause  a  plane  wave,  incident  from  below,  to  produce  interference 
patterns  at  the  surface.  Similar  effects  may  explain  the  variation  in 
character  of  the  seismic  arrivals  at  individual  subarrays. 

The  identification  of  reflections  from  depths  corresponding  to 
the  Mohorovi^id  discontinuity  has  an  important  bearing  on  the  controversy 
surrounding  the  nature  of  the  crust-mantle  boundary.  Almost  since  the 
time  of  A.  Mohorovicid' s  discovery  in  1909  of  a  refracting  interface  at 
depths  of  a  few  tens  of  kilometers,  two  hypotheses  have  been  advanced 
concerning  the  mineralogy  at  the  base  of  the  crust.  One  holds  that  the 
Moho  is  a  chemical  boundary  separating  the  basic  rocks  (predominantly 
basalt)  of  the  lower  crust  from  the  ultrabasic  rocks  (peridotite)  of  the 
upper  mantle.  The  second  assumes  that  the  boundary  represents  an  iso¬ 
chemical  phase  change  from  basalt  to  eclogite.  In  recent  years,  Lovering 
(1958),  Kennedy  (1959),  MacDonald  and  Ness  (1960),  Stishov  (1963)  and 
others  have  advocated  the  hypothesis  of  a  phase  change  for  the  contin¬ 
ental  crust-mantle  boundary.  Bullard  and  Griggs  (1961),  Nakamura  and 
Howell  (1964) ,  Ringwood  and  Green  (1964) ,  Ringwood  (1966)  and  others 
have  opposed  this  idea  and  argue  in  favor  of  a  chemical  transition.  If 
the  Moho  is  a  phase  transformation,  it  is  generally  accepted  that  the 
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velocity  transition  would  occur  over  a  vertical  distance  of  at  least 
a  few  kilometers  and  possibly  as  much  as  15  kilometers.  This  is  contrary 
to  the  evidence  presented  in  this  study  in  which  near- vertical  incidence 
reflections  are  obtained  from  depths  corresponding  to  that  of  the  re¬ 
fracting  discontinuity.  Such  data  support  the  hypothesis  that  the  base 
of  the  crust  is  a  sharp  boundary,  less  than  one  kilometer  thick,  repre¬ 
senting  a  change  in  chemical  composition. 

This  study  contains  the  first  attempt  at  obtaining  the  specific 
attenuation  factor,  Q,  as  a  function  of  depth  from  reflected  seismic 
waves.  On  the  basis  of  a  comparison  of  autopower  spectral  densities  of 
field  records  and  synthetic  seismograms,  an  acceptable  Q  structure 
appears  to  require  considerable  variation  in  the  sedimentary  layers,  but 
having  an  average  value  of  nearly  300.  In  the  crystalline  basement, 
the  magnitude  of  Q  is  approximately  1500  and  probably  increases  with 
depth.  In  order  to  provide  values  which  have  more  reliability  and 
precision,  wide-band  recording  of  reflections  from  within  the  sediments, 
with  similar  instrumentation  and  field  procedures  as  that  used  for  the 
deep  reflections,  would  be  valuable.  Indeed,  it  is  likely  that  one  may 
be  able  to  determine  the  distribution  of  Q  with  depth  and  also  map 
lateral  variations  from  an  average  model. 

Initial  attempts  at  determining  the  nature  of  deep  crustal  reflec¬ 
tors  has  met  with  more  success  than  was  originally  anticipated.  Auto¬ 
power  spectra  and  synthetic  seismograms  have  been  utilized  to  specify 
models  of  transition  zones  with  acceptable  characteristics.  The  para¬ 
meters  used  for  the  study  included  the  shapes  of  the  spectra,  the 
relative  rates  of  attenuation  of  frequency  components,  the  forms  of 
reflected  wavelets  for  a  pulse  input  and  the  ratios  of  the  power  in  the 
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spectra  of  shallow  reflections  to  that  in  the  deeper  ones.  Models  incor¬ 
porating  a  first-order  discontinuity  in  velocity,  a  linear  increase  of 
velocity  with  depth  or  step  increases  in  velocity  did  not  satisfy  all 
the  observational  data.  Thin  layers,  less  than  200  meters  thick,  of 
alternating  velocity  increases  and  decreases  over  a  total  depth  extent 
of  less  than  one  kilometer  provided  a  model  which  was  acceptable  on 
the  basis  of  observed  results. 

Thus  the  import  of  seismic  reflection  crustal  studies  is  already 
being  realized.  The  method  can  provide  the  details  of  deep  crustal 
structure  necessary  for  the  interpretation  of  a  much  wider  range  of  data. 
It  is  hoped  that  this  research  has  contributed  significantly  toward  the 
establishment  of  the  seismic  reflection  technique  as  a  recognized  and 
viable  method  of  crustal  research. 
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APPENDIX 

FORTRAN  SOURCE  LISTINGS  OF  COMPUTER  PROGRAMS 

The  velocity  filter  program 

The  velocity  filter  program  is  designed  to  process  deep  reflection 
seismograms  with  a  fan-pass  filter.  The  logic  of  the  program  follows  the 
algorithm  developed  by  Treitel  et  al .  (1967).  The  input  data,  comprised 
of  amplitude  values  from  two  seismograms,  are  processed  using  three  dif¬ 
ferent  filter  characteristics.  This  data  and  the  three  sets  of  filtered 
output  data  are  plotted  by  the  use  of  CALCOMP  plotter  subroutines.  As 
well,  the  output  data  are  written  on  tape. 

Execution  time  for  the  program  is  short.  If  twelve  vectors  of 
960  points  are  processed  by  three  different  filters  which  use  eight  of 
these  data  sets  per  single  output  vector,  each  filter  thus  producing  a 
total  of  five  output  traces,  the  execution  time  is  0.92  minutes  on  the 
IBM  360/67  computing  system. 
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C 

C  THIS  PROGRAM  VELOCITY  FILTERS  SEISMIC  DATA  ACCORDING  TO  THE 
C  ALGORITHM  DEVELOPED  BY  TREITEL  FT  AL .  (GEOPHYSICS,  V,  32, 

C  P.  739-800,1967}. 

C  THE  BLOCK  LENGTH  OF  THE  INPUT  DATA  IS  1200  WORDS/BLOCK 
C  NOTE. .... .FOR  VELOCITY  FILTERED  DATA  ONLY  THE  OUTPUT  TAPE  COMPRISES 

C  DATA  FOR  WHICH  THE  BLOCK  LENGTH  IS  200  WORDS/BLOCK 

C 

1000  FORMAT  (1615) 

1001  FORMAT  (8F10.5) 

1002  FORMAT  <»0  THE  TOTAL  NUMBER  OF  PLOTS  FOR  INCLUDING  ON  PLOT  REQUEST 

1  SLIP  IS", 15) 

1003  FORMAT  tlHO, ’POINTS  1  TO®, 15,®  FOR  THE  DATA  POINTS®, IS, 1  TO®, 15,® 
IGF  INPUT  TRACE  ',12,®  WHICH  ARE  TO  BE  FILTERED') 

1004  FORMAT  (1H0,  10 { IX , E  1 1 . 4 J  ) 

1005  FORMAT  { 1 H  0 , 'POINTS  1  TO', 15,®  FOR  THE  ® ,12,'TH  VELOCITY  FILTERED 
10UTPUT  TRACE  WHOSE  POSITION  IS  * ,F4.1,®  RELATIVE  TO  POSITION  I*) 

1006  FORMAT  (1H0,I4,®  BLOCKS  OF  200  WORDS  PER  BLOCK  WERE  COPIED  ONTO  TH 
IE  NEW  TAPE5 ) 

1007  FORMAT  (  1H1 , • BLOCKS  *,I5,®  TO  ®, I  5,® (200  WORDS  PER  BLOCK)  CONTAIN 
1 THE  * ,12,'TH  VELOCITY  FILTERED  OUTPUT  TRACE  FROM  THIS  DATA') 

100  8  FORM  AT  ( 9 1  5 , 5 A4 , 3F 5 . 1 ) 

1009  FORMAT  (1H0,5A4,®  HAVE  BEEN  VELOCITY  FILTERED  TO  PROV IDE 1 , 1 3 , '  TRA 
ICES  OF',  15,®  POINTS/TRACE,  THE  LAST  DATA  ON  TAPE*//®  FILTER  CENTER 
2ED  AT  ®,I2,®  DELTA  T/TRACE  WITH  A  CUTOFF  OF  +  OR  -  ',12,®  DELTA  T/ 
3TRACE® ) 

G 

1010  FORMAT  {  1H 1 »  *  NUMBER  OF  INPUT  TRACES  EQUALS  1  , I 3//2X , * NO.  OF  TRACE 
IS  USED  IN  FAN  PASS  OPERATION  EQUALS  1  ,  I2//2X, 'PASS  BAND  CENTERED  A 
2B0UT  *  »  I  2  ,  *  DELTA  T/TRACE.  CUTOFF  APPARENT  VELOCITY  IS  ®,I2,®  DEL 
3T  A  T/TRACE* //2X , 1  NO .  CF  POINTS  TO  BE  FILTERED  EQUALS  * ,  I  5//2X , '  ST  A 
4RTING  BLOCKS  ON  TAPE  FOR  INPUTTING  DATA  FROM  2  RECORDS  ARE  *,15,® 
SAND  »  ,  I  5/  f2X  ,  *FCR  DATA  READ  INTO  MATRIX  XDATA,  THE  STARTING  POINT 

6 FOR  VELOC  TY  FILTERING  IS®, 15) 

C 

COMMON/ST  )R  AG/ XD  AT  A ( 12 » 2400 ) , BLKST ( 2 } 

COMMON/ST  IRE/XFILT (5,1600),  SEPT,  T  S  T  ART ,  TINCH,  SCALE,  COPS, 

1  RECN015)  ,  CENTER,  NT AU 

DIMENSION  SUM { 1 600 ) ,  Y(16C0),  Z11600),  SPARE(200),  BUFF ( 1536} 
INTEGER  B  .KST,  PTSTAR,  CENTER,  START,  ROWST,  PTEND,  CENT ( 3 )  , T AU ( 3  ) 
EQUIVALEN  IE  I SUM( 1 ) , SPARE ( 1 ) ) 

XPOS  =  2  .  ) 

YPOS  =  25  >0 
IP  =  0 
NREAD  =  0 

READ  (5,  1000)  M,  NBLK ,  NT  R ,  NOUT,  NCOPY*  NPRIMT,  JPLOT 

M  MUST  BE  EVEN  (8,  10  OR  12)  AL  0  GIVES  THE  NO.  OF  TRACE*  TO  BE 

USED  IN  THE  FAN  PASS  OPERATION 

NBLK  NO.  IF  BLOCKS  OF  DATA  OF  6  X  200  POINTS  IN  EACH  RECORC 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NTR  NUMBER  OF  INPUT  TRACES  WHICH  ARE  TO  BE  FILTERED 

NOUT  THE  NUMBER  OF  OUTPUT  TRACES  EXPECTED 

NCOPY  NO.  OF  BLOCKS  OF  200  WORDS/BLOCK  TO  COPY  FROM  UNIT  3  ONTO  2 

NPRINT  NO.  OF  DATA  VALUES  TO  LIST 
JPLO.T  =  0  AMD  NO  PLOTTING  IS  DONE 

=  1  AND  THE  INPUT  AND  VELOCITY  FILTERED  TRACES  ARE  P L 0 T T E 0 

READ  (5,10011  SEPT,  TINCH,  COPS,  SCALE 

THIS  READS  IN  PARAMETERS  REQUIRED  FOR  THE  PLOTVF  ROUTINE 
SEPT  SEPARATION  IN  INCHES  BETWEEN  TRACES 
TINCH  NUMBER  OF  INCHES  OF  PAPER  PER  SECOND 

COPS  NUMBER  OF  CONVERSIONS  PER  SECOND 

SCALE  SCALE  FACTOR  TO  REDUCE  MAGNITUDE  FOR  PLOTTING  IN  INCHES 


IF  (JPLOT.EQ.  1)  CALL  PLOTS  ( BUFF ( 1 } , 6 144 3 
NTR  PT  =  0 

IF  (NCOPY. EQ.O)  GO  TO  4 
DO  2  J-l  , NCOPY 
READ  ( 3 , ERR=3 )  SPARE 
GG  TO  2 

3  READ  (3)  SPARE 
2  WRITE  (2)  SPARE 

4  WRITE  (6,10061  NCOPY 
NWRIT  =  NCOPY 

1  READ  (  5,  1008  )  ISTOP,  BLKST(l),  BLKST ( 2 ) ,  (CENT (II, 1  =  1, 3),  (  TALK  I  )  , 
11  =  1,3),  (RECNO(  I  ) ,  1  =  1  ,5) ,  RST  ART ,  TST  ART •  TEND 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ISTOP 
BLKST 
CENT ( I  1 


T  AU ( I  1 


RECNO 
R  START 
T  START 
TEND 
TST  ART 


=  099  IN  COLUMNS  3-5  STOPS  THE  RUN 
STARTING  BLOCK  FOR  EACH  GROUP  OF  6  TRACES  WHICH  MAY  BE  USED 
IS  THE  MOVEOUT  IN  SAMPLING  INCREMENTS/ TRACE  ABOUT  WHICH  THE 
PAS'S  BAND  IS  CENTERED.  FOR  2  APPARENT  VELOCITIES,  CENTER 
IS  DETERMINED  FROM  { V2  +  VI 1/2. 

IS  THE  MOVEOUT  IN  SAMPLING  INCREMENTS/TRACE  WHICH 
APPROXIMATES  THE  CUTOFF  APPARENT  VELOCITY.  FOR  2  VELOTITIES, 
NTAU  IS  DETERMINED  FRGM  ABSCV2'-  Vll/2. 

ALPHAMERIC  IDENTIFICATION  FOR  DATA  BEING  FILTERED 
TIM!  CORRESPONDING  TO  POINT  I  OF  THE  RECORDS 


TIM! 

TIM) 


AT  WHICH  VELOCITY  FILTERING  AND  PLOTTING  BEGIN 
AT  WHICH  VELOCITY  FILTERING  AND  PLOTTING  END 


-  TEND  MUST  BE  LESS  THAN  OR  EQUAL  TO  13.2  SECONDS 


IF  (ISTOP  EQ.9195  GO  TO  110 

PT  =  (TST..RT  -  R  ST  ART )  *  COPS 

IF  ((PT  -  IFIX(PT) l.GE.0.5)  PT  =  Pr  +  1.0 

PTSTAR  =  IFIX(PT)  +  1 

PT  =  (TEND  -  RST ART  1  *  COPS 

IF  ((PT  -  IFIX(PT)  l.GE.0.5)  PT  =  P T  +  1.0 

NPEND  =  I!  I  X  (  PT  ) 

NPTS  =  NP  iND  -  PTSTAR  4-  1 
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CALL  READ3 ( NBLK , NREAD ) 

C  THIS  FEEDS  INTO  XDATA  THE  TRACES  TO  BE  VELOCITY  FILTERED 
C 

NP  =  PT STAR  +  NPRINT  -  1 
DO  5  J- 1 ,  NT  R 

WRITE  (6,1003)  NPRINT,  PTSTAR,  NP END ,  J 
5  WRITE  (6,1004)  ( XDATA ( J , J A) , JA=PT STAR ,NPJ 

IF  ( J PLO T  •  EG . 1 )  CALL  PLCTVF  1 NTR, NPTS , 1 , NOUT , PTSTAR ,  0,NTR  PT, 
1XP0S,  Y  P  C  S ,  IP) 

THIS  USES  THE  CALCCMP  PLOTTER  TO  PLOT  THE  INPUT  TRACES  TO  THE  V.F. 

DO  100  K F  =  1,3 
CENTER  =  CENT (  KF ) 

NT  AU  =  TAU(KF) 

WRITE  (6,1010)  NTR,  M , CENTER, NTAU, NPTS , BLKST ( 1 ), BLKST ( 2 ), PTST AR 
L  =  M-l 

MC-  =  L* I ABS( CENTER) 

NEND  =  NPEND  -  (L+l)*NTAU/2 
START  =  PTSTAR  +  ( L~ 1 ) *NT AU /2 
IF  ( CENT  ER . LT . 0 )  START  =  START  -5-  MC 
IF  (CENTER.  GT  .0 )  NEND  =  NEND  -  MC 
C  THE  STARTING  AND  ENDING  TIMES  HAVE  BEEN  CHANGED  SO  THAT  ONLY  POINTS 
C  IN  THE  TIME  INTERVAL  T START  TO  TEND  WILL  BE  USED  WHEN  THE  DATA 
C  POINTS  ARE  TIME  SHIFTED. 

C 

A  =  0.65465 
B  =  0.93632 
C  =  0.13091 
ROWST  =  0 

20  ROWST  =  ROWST  +  1 
DO  10  1=1, 1600 
10  SUM {  I)  =  0.0 

DO  50  N=  ST  ART ,  NEND 

NSN  =  N  -  PTSTAR  +  1 

IROW  =  RCl  ST 

DO  40  J  =  1  <  M 

JA  =  J-l 

MU  =  2  * J  A  -  L 

MUDE  =  (  M'-l  )*NTAU/2 

MUAD  =  ( Ml'+l  )*NTAU/2 

JB  =  CENT!  R*JA 

I  COL 1  =  N  -  MUDE  +  JB 

IC0L2  =  N  +  MUAD  +  JB 

<  DIFF  =  ( XI  ATA( IROW , IC0L1 )  -  XDATA (  [ROW , IC0L2 )) /MU 

SUM(NSN)  =  SUM(NSN)  *  DIFF 
40  IROW  =  IROW  +  1 
50  CONTINUE 

NS T 1  =  ST)  RT  -  PTSTAR  +  1 
NEN 1  =  N  El  D  -  PTSTAR  +  1 

C  THE  VECTOR  »IUM*  HAS  NOW  BEEN  FILLED  IITH  THE  WEIGHTED  SUM 
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C  OF  THF  TIME  SHIFTED  TRACES, 
r 

DO  55  1=1,1600 
Y(I)  =  0.0 
Z  (  I  )  =  0.0 

55  XFILKROVtST,  I)  =0.0 
DO  60  N  =  N  S  T 1 , N  E  N  I 
NA  =  N-l 
NO  =  N-2 

Y ( N }  =  SEMINA)  -  A  *  S  U  M  (  N  8  )  +  B  * Y  C  N  A )  -  C*Y(NB) 

NN  =  NEN 1  -  K  +  NST 1 
NC  =  NN+1 
ND  =  N  N  + 2 

60  Z(NN)  =  SUM  INN)  -  A*SUM(NC)  +  8*Z(NC)  ~  C*Z(ND) 

DO  70  N=NST  1  ,NEN 1 

70  XE  ILT  (  ROViST  ,N  )  =  0 .2026424*  (  Z  {  N  )  -  Y  {  N )  ) 

LL  =  NEN  1  4  ? 00  -  MOD [ NEN 1 , 200 J 

C  THE  RECURSIVE  FILTER  OPERATOR  HAS  JUST  BEEN  APPLIED  TO  THE  DATA 
C  X  F  ILT { RGWST ,  )  CONTAINS  THE  (ROWST)TH  VELOCITY  FILTERED  OUTPUT  TRACE 

C 

IF  (ROWST.LT.NGUT)  GO  TO  20 

IF  (JPLOT.EQ. 1 )  CALL  PLOTVF  ( ROW S T , NP TS , 2 » NOUT , 1 , Kf ,  NT R  PT ,  XPOS, 
1  YPGS,  IP) 

FM  =  FLOAT (M/2) 

DO  80  J=1,N0UT 

POS  =  FM  4-  0.5  4"  FLOAT  {  J—  1 ) 

LL  L  =  LL/200 

NW  =  NWRIT  4-  1 

NWRIT  =  NWRIT  4-  LL  L 

WRITE  (6,3007)  NW,  NWRIT,  J 

WRITE  (  6,1005)  N  PR  I  NT ,  J,  POS 

DO  75  JA  =  1 ,  LLL 

JC  =  200*.  A 

J3  =  JC  -  199 

75  WRITE  (2)  ( XF  ILK J , JD ) , JD=JB» JC) 

80  WRITE  (6,  004)  { XF ILT (J , J A)  , J A= 1 , NPR I  NT  ) 

WRITE  (6,  .009)  R  EC  NO ,  NOUT,  LL,  CENTER,  NT  AU 
100  CONTINUE 
GO  TO  I 

110  IF  (JPLOT.EQ. 0)  GO  TO  120 

CALL  PLOT  (XPOS, YPOS-28. 0,-3) 

IP  =  IP  4-  2 

,  CALL  PLOT  (0.0,  0.0,  999) 

WRITE  (6,  002)  IP 
120  END  FILE  2 
REWIND  2 
STOP 
END 


o  o  o 
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SUBROUTINE  READ3  ( NBLK , NRE AD ) 

C  THIS  SUBROUTINE  READS  12  CHANNELS »  6  FROM  EACH  OF  2  DIFFERENT 
C  RECORDS  ?  CF  UP  TO  2400  POINTS  INTO  XDATA 
C 

C  MOUNT  DATA  TAPE  ON  FTC04 

C  NBLK  IS  THE  NO.  CF  BLOCKS  OF  6  X  200  POINTS  TO  BE  READ  IN 
C  BLKST  IS  THE  STARTING  BLOCK  FOR  EACH  GROUP  OF  6  TRACES  WHICH  ARE  USED 

C 

COMMGN/STORAG/XDATAI 12 , 2400) ,  BLKST { 2 ) 

INTEGER  BLKST 
NPTS  =  200*NBLK 
DO  30  K A  =  1 ,2 

NSKIP  =  BLKST [ K A )  -  NREAD  -  1 
CALL  SKIP  (NSKIP,  4) 

NREAD  =  NREAD  4-  NSKIP 
10  JB  =  6*KA 
JA  =  JB  -  5 
DO  20  J=l, NPTS, 200 
JC  =  J  +  199 

READ  I  4  i  £RR=15 )  ( ( XDATA ( J D, JE ) t JD= JA, JB) , JE= J T JC ) 

GO  TO  20 

15  READ  (4)  ( (XOATA(JO, JE) tJD=JA, JB)  ,JE=Jt JC) 

20  CONTINUE 

NREAD  =  NREAD  +  NBLK 
30  CONTINUE 
RETURN 
END 


SUBROUTINE  SKIP  ( I SK  I  P , I T APE ) 

THIS  SUBROUTINE  IS  USED  TO  SKIP  DOWN  A  TAPE  TO  THE  DESIRED  RECORD 

IF  (ISKIP.EQ.O)  RETURN 
IF  ( ISKIPoLT .0)  GO  TO  1 
DO  2  II  =  1  ,  ISKIP 
READ  (ITPE,  ERR  =  5) 

GO  TO  2 

5  READ  (  IT  A  3E ) 

2  CONTINUE 
RETURN 

1  ISKIP  =  -  ISKIP 
DO  3  1  =  1  ,  ISKIP 
BACKSPACE  ITAPE 

3  CONTINUE 
RETURN 
END 


o  n  n  r> 
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SUBROUTINE  PLOTVP  ( NTP , NPTS , I NDEX , NOUT , NPST AR , KF , NTR  PT , XPOS, YPOS , 

1  IP) 

THIS  SUBROUTINE  PLOTS  THE  INPUT  DATA  FOR  THE  VELOCITY  FILTER  AND  ALL 
DATA  WHICH  HAS  BEEN  FILTERED 

CGMMON/STOR AG/XDAT A( 12,2400) ?  8  L  K  S  T ( 2 ) 

COMMON/STORE/XFILT 15,  1600),  SEPT,  T  ST  ART ,  TINCH,  SCALE,  COPS, 

1  RECNO ( 5  ) ,  CENTER,  NTAU 
INTEGER  CENTER 

WRITE  (6,1001)  N PT S, NTP, SEPT, TINCH, COPS, TST ART, SCALE 
1001  FORMAT  I  1  HO , 1  FOR  PLOTTING  THE*,  15, 5  POINTS  FOR*, 13,*  TRACES,  THE  T 
1RACE  SEPARATION  I S * , F 5  * 2 / /4 X ,  s THE  PAPER  SPEED  IS  8,F4«2,*  IN. /SEC 
2 AN  D  THE  DIGITIZING  RATE  IS  * ,F7.2,e  CONV/SEC * // 4X, 5  THE  BEGINNING  T 
3  I M E  IS  *  , F5 . 2 , *  AND  THE  SCALING  FACTOR  IS  S,F7.1) 

NP END  -  NPSTAR  +  NPTS  -  1 
T  =  TINCH/COPS 
XLEN  =  FLOAT  (NPTS)  *T 
A8SINT  =  l.O/TINCH 
MNTR  =  MCD  (NTR  PT ,  2) 

GO  TO  (20,30),  INDEX 
20  DO  27  J=1,NTP 

CO  22  JA  =  NPSTAR, NPENO 
22  XDAT  A ( J , JA )  =  XDAT A ( J , JA ) /SCALE 
27  CONTINUE 

CALL  PLOT  ( XPOS, YPOS, -3) 

IP  =  I  P+1 

CALL  PLOT  (0. 0,-1. 0,-3) 

IP  =  IP+I 

CALL  AXIS  (0.0, 0.0,*  TIME  IN  SECONDS  * , +1 5 , XLEN , 0 . 0, TSTART , AB S I  NT , 

110.0) 

YD  =  3.5  -f  5. 5  *  SEPT 

CALL  SYMBC  L  (- 1 . 0 , -YD , 0 . 2  , RECNO , 90 . 0 , 20 ) 

YDADD  =  Yl  +  0.50 

CALL  SYMB(  L  (-0.5, -YDADD  ,0 . 15  ,  *  D  I  STANCE  BETWEEN  TRACES  IS  .29  KM. 
I  * ,90.0,34' 

GO  TO  39 

30  DO  10  J=1  NTP 

DO  5  JA  =  NPSTAR, NPENC 
5  XFILTIJ.J/ )  =  XF ILT( J ,JA) /SCALE 
10  CONTINUE 

IF  (MNTR. HE. 0)  GO  TO  36 
’  33  CALL  PLOT  ( -XL EN  , 0 . 0 , -3 ) 

IP  =  IP  +  1 

36  IF  (KF.NE.l)  GO  TO  39 
CALL  PLOT  (0.0, -2. 0,-3) 

IP  =  IP  +  1 

CALL  AXIS  ( 0.0, 0.0, *  TIME  IN  SECOND S * ,  +  1 5 , XL EN , 0 . 0 , T ST  ART , A B 5 1  NT , 

110.0) 


. 
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39 


40 

41 


42 

43 


44 


45 

46 

47 


CALL  SYMBOL  (-1  .0,-0. 7, .20, ’VELOCITY®  ,90.0,8) 
CALL  SYMBOL  ("0.7, -0.7,0,20?  *  FILTERED*  ,90.0,8) 
CALL  SYMBOL  { -0 . 3 , -0 , 4 , 0 . 20 , *  DAT A  * , 90  .  0 , 4 ) 

IF  (KF.EQ.O)  GO  TO  44 

YD  =  2.1  +  FLOAT { NGUT-1 ) *SEPT 

FNUM1  -  FLOAT  (CENTER  -  NTAUJ/COPS 

IF  ( FNUM  l.NE.O.  )  GO  TO  40 

VEL  I  =  1  .0 

GO  TO  41 

VEL 1  =  0.293/FNUM1 
FNUM2  =  FLOAT  ( CENTER+NT  AU) /COPS 
IF  (FNUM2.NE.0. )  GO  TO  42 
VEL2  =  1  .0 
GO  TO  43 
VEL2 
CALL 
CALL 
CALL 
YP  = 

CALL 
CALL 
YP  = 

CALL 
CALL 
YP  = 

CALL 
CALL 
CALL 
CALL 
YP  = 

CALL 
CALL 
YP  = 

CALL 
CALL 
YP  = 

CALL 

DO  80  J=:  y  NT  P 
MJ  =  MOD  J  y  2 ) 

IF  ( J.NE. 1)  GO  TO  45 
CALL  PLO'  ( 0  .0,-2. 0,-3) 

IP  =  IP  -  1 

GO  TO  46 

IF  (MJ.EQ.O)  GO  TO  60 
CALL  PLOT  (— XLEfSy-SEPT  >  —  3  ) 

IP  =  IP+1 

CALL  SYMBOL  (0. 0,0.0, .15,04,0.0,- ♦) 

GO  TO  (47,49),  INDEX 
CALL  PLOT  (0.0, XDATA(J,NP STAR) , 3 ) 

DO  48  JA  =  NPST AR, NPEND 


=  0.293/FNUM2 

SYMBOL  (-1.0,- YD,  .  15,  'VI1  ,90.0,2) 
SYM80L(-0.9,-YD+.25, . 1 0 , 1 APP * , 90 . 0 , 3 ) 
WHERE  { XP, YP, FCTR) 

YP  4  0.2 

SYMBOL  (-1.0, YP,  .15,  '  IS  ',90.0,4) 
WHERE  ( XP , Y  P , FC  TR ) 

YP  4  0.2 

NUMBER  (-1 .0,YP,  .15, VEL l, 90.0, 1) 

WHERE  ( XP , YP , FCTR) 

YP  4  0.2 

SYMBOL  (-1.0, YP,  .15,  'K/S*  ,90.0,3) 
SYMBOL (-0.5, -YD, .15, 8V2',9Q.0,2) 
SYMBOL (-0.4, -YD 4.2 5, .10, *APP» ,90.0,3) 
WHERE  ( XP , Y  P , FCTR) 

YP  4  0.2 

SYMBOL  (-0.5, YP,  .15, 8  IS  *,90.0,4) 
WHERE  (XP,YP,ECTR) 

YP  h  0.2 

NUMf  ER  (-0.5,YP,  .  15 , VEL2 , 90 .0 , 1) 

WHEJ E  ( XP, YP, FCTR) 

YP  -  0.2 

SYMCOL  (-0.5,YP,  .15,  'K/S* ,90,0,3) 


TSCALE  =  T  *  FLOAT ( JA-NPSTAR+1 ) 

48  CALL  PLOT  ( TSCALE * XDATA ( J , J A ) , 2 ) 

GO  TO  75 

49  CALL  PLOT  ( 0. 0 , XF I LT (J * NPST AR > , 3 ) 

DO  50  JA  =  NPSTAR, NPEND 

TSCALE  =  T  *  FLOAT (JA-NPSTAR+1) 

50  CALL  PLOT  ( T SC ALE , XF I LT ( J  ,  J A ) , 2  ) 

GO  TO  75 

60  CALL  PLOT  ( XLEN  ,-SEPT  f-3 ) 

IP  =  IP+1 

CALL  SYMBOL  ( 0  .  C , 0  . 0 ,  .  1 5 , 04 , 0 . 0 ,-4 ) 

GO  TO  161,64) »  INDEX 

61  CALL  PLOT  ( 0*0, X DATA! J » NP END )  , 3 ) 

DO  62  JA  =  NPSTAR,NPEND 

TSCALE  =  -  T  *  FLOAT! JA-NPSTAR+1 ) 

JCC  =  NP  END  -  JA  +  NP  ST  AR  -  1 
IF  ( JCC • LT .NPST AR )  JCC  =  NPSTAR 

62  CALL  PLOT  (TSCALE,  XD ATA { J , JCC ) , 2 ) 

GO  TO  75 

64  CALL  PLOT  (0*0, X  FILE ( J , N  P  EN  D ) ,  3  ) 

DO  65  JA  =  NPSTAR,NPEND 

TSCALE  =  -  T  *  FLOAT! JA-NPSTAR+1) 

JCC  =  NP  END  -  JA  +  NPSTAR  -  1 
IF  ( JCC. LT. NPSTAR)  JCC  =  NPSTAR 

65  CALL  PLOT  (TSCALE,  XF JLT ( J , JCC )  , 2 ) 

75  CALL  PLOT  ( TSCALE, 0.0,3) 

CALL  SYMBOL  ( TSCAL E, 0 . 0 , .  15 , 04, 0 . 0 ,-4 ) 

80  CONTINUE 

GO  TO  (31 ,84) ,  INDEX 

81  DO  83  J=  1  ?  NT  P 

DO  82  JA  =  NPST  AR , NP  END 

82  XD AT  A  {  J , JA )  =  XDAT A ( J , JA) *SCALE 

83  CONTINUE 
GO  TO  92 

84  DO  90  J=]  , NT P 

DO  85  JA  =  NPSTAR, NPEND 

85  XFILT(J,v  A)  =  XF ILT ( J , JA) TSCALE 
90  CONTINUE 

92  ’WRITE  (6,1000)  IP 

1000  FORMAT  (1  HO,* NUMBER  OF  TIMES  ORIGIN  IS  SET  TO  ZERO  IS 
IF  ( KF.NI  .3 )  GO  TO  100 
IF  (MNTR.EQ.O)  GO  TO  95 
XPGS  =  XI EN  +  5.0 
GO  TO  97 
95  XPOS  =5.0 

97  YPOS  =11.0  +  FLOAT (3*NTP  +  8)  *  SEPT 
ICO  NTR  PT  =  NT P 
RETURN 
END 
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The  synthetic  seismogram  program 

The  principal  basis  of  the  synthetic  seismogram  program  is  an 
iterative  formula  for  calculation  of  the  transfer  function  of  a  layered 
medium  in  which  the  velocity  may  vary  linearly  with  depth  in  any  layer. 
The  formula  is  from  a  development  by  Berryman  et  al.  (1958)  which  has 
been  rederived  in  order  that  the  effects  of  attenuation  may  be  included. 
The  computed  transfer  function  is  plotted. 

A  choice  of  two  different  input  waveforms  is  provided  in  the 
program.  One  of  these  is  the  form  of  Gram-Charlier  wavelet  used  by 
Berryman  et  al .  (loc.  cit.).  By  proper  specification  of  the  input 
parameters,  the  principal  frequency  component  of  the  wavelet  can  be 
controlled.  The  second  form  of  input  pulse  is  specified  by  its 
amplitude  versus  frequency  characteristics,  assuming  a  constant  phase 
of  zero.  An  inverse  Fourier  transform  of  this  is  taken  to  produce  an 
approximation  to  a  delta  function  centered  at  zero  time.  The  pulse 
is  then  delayed  and  a  Fourier  transform  of  the  delayed  pulse  is  computed 
to  obtain  both  real  and  imaginary  parts  of  the  pulse  spectrum. 

An  inverse  Fourier  transform  of  the  product  of  the  complex 
pulse  spectrum  and  the  transfer  function  yields  the  amplitude  values 
of  the  synthetic  seismogram.  These  are  plotted  on  a  specified  time 
scale. 

A  facility  is  included  to  allow  computation  of  the  autopower 
spectra  from  up  to  four  different  time  intervals  on  the  synthetic 
seismogram.  The  autopower  spectral  density  of  the  first  interval  is 
plotted  over  the  entire  range  from  zero  to  the  Nyquist  frequency.  Any 
other  spectra  computed  are  plotted  on  one  additional  graph  and  for  the 
frequency  range  from  5  to  25  Hz. 
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The  execution  time  of  the  program  is  controlled  by  the  calculation 
of  the  transfer  function  of  the  layered  medium.  A  large  number  of  layers 
coupled  with  the  necessity  of  having  small  frequency  increments  causes 
the  transfer  function  to  be  computationally  expensive.  The  generation 
of  the  synthetic  seismograms  used  in  this  thesis  required  considerable 
computing  time.  Because  a  velocity  log,  averaged  over  20  foot  intervals, 
was  used  to  model  the  sedimentary  layers,  about  242  layers  were  intro¬ 
duced  for  this  part  only.  Depending  on  the  complexity  of  the  remaining 
part  of  the  specified  layered  medium,  approximately  250  layers  would 
have  to  be  considered.  In  order  to  generate  a  seismogram  at  least  18 
seconds  long  required  a  small  frequency  increment  since  the  inverse  of 
this  increment  is  the  maximum  time  interval  which  can  be  calculated 
without  aliasing  problems.  The  frequency  increment  is  constrained  by 
the  specified  Nyquist  frequency  and  the  requirements  of  the  fast  Fourier 
subroutines  which  use  a  number  of  points  that  must  be  a  power  of  two. 

For  most  calculations,  the  writer  used  a  Nyquist  frequency  of  80  and 
N/2  =  11  resulting  in  frequency  increments  of  80/211.  On  the  basis  of 
the  quantities  mentioned,  the  total  execution  time  of  the  synthetic 
seismogram  program  on  the  IBM  360/67  computing  system  was  about  16 
minutes.  On  the  other  hand,  for  a  ten- layer  case,  a  Nyquist  frequency 
of  80  and  N/2  =  10,  the  execution  time  was  less  than  1  minute. 
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THIS  PROGRAM  IS  DESIGNED  TO  PRODUCE  A  SYNTHETIC  SEISMOGRAM  AT 


VERTICAL  INCIDENCE  ALLOWING  FOR 
THE  PROCESS,  IT  PLOTS  THE  INPUT 


ATTENUATION  WITH  DEPTH,  IN 
PULSE  y  ITS  FOURIER  TRANSFORM, 


THE  REFLECTION  COEFFICIENT  OF  THE  HORIZONTALLY  LAYERED  MEDIA, 


THE  FINAL  SYNTHETIC  SEISMOGRAM  AND  AUTOPOWER  SPECTRA  OF  SPECIFIED 
TIME  INTERVALS  OF  THE  SEISMOGRAM, 


C 


999  FORMAT  (16F5.0) 

1000  FORMAT  (4I5,12F5.0) 

1001  FORMAT  {’1  THE  SYNTHETIC  SEISMOGRAM  AMPLITUDE  VALUES1//) 

1002  FORMAT  { 10 ( IX , E 1 1 . 3  )  ) 

1003  FORMAT  (*1  THE  NUMBER  OF  POINTS  IN  THE  TIME  SERIES  IS  *,15// 

1*  THE  NUMBER  OF  POINTS  BETWEEN  FREQUENCIES  OF  0  AND  NYQUIST  IS 
2X5//S  THE  NYQUIST  FREQUENCY  FOR  THIS  RUN  OF  THE  PROGRAM  IS  8,F7.1 
3// *  THE  CORRESPONDING  DIGITIZING  INTERVAL  IS  * , F7.3// 

4*  THE  APPROXIMATE  PERIOD  OF  THE  INPUT  PULSE  IS  * ?F7.3// 

5*  FOR  PLOTTING  TIME  SERIES,  THE  NO,  OF  INCHES/SECOND  IS  *,F7.3// 
6*  THE  AMOUNT  OF  TIME  TO  BE  PLOTTED  FOR  THF  S.  $.  IS  NF7.3) 

1004  FORMAT  ( 2 [ 5 X , F 1 0 . 3 , 5X  ,  E 1 5 , 5 ) ) 

1005  FORMAT  {'1  FREQUENCY  AND  AMPLITUDE  VALUES  FOR  THE  A.P.  SPECTRA  OF 
1THE  S.  S.  FROM  f,F5.2,‘  SEC.  TO  • , F5.2,f  SEC.') 

1006  FORMAT  (H  THE  GAIN  CONTROL  FACTOR  AS  A  FUNCTION  OF  TIME  FOLLOWS.  * 
l//s  TIME  STARTS  AT  0,  WITH  AN  INCREMENT  BETWEEN  VALUES  OF  TWICE  * 

‘  2  ,  F6 , 3  , 6  SEC.*//) 

1007  FORMAT  1*0  INCIDENT  PULSE  SPECTRUM  IS  ZERO  FROM  0  HZ.  TO  *  ,F5.1, 
1*  REACHES  UNITY  AT  *,F5.1,*  AND  SYMMETRICALLY  RETURNS  TO  ZERO  AT  * 
2fF5*l//'  THE  GAIN  IS  (,F8.1,1  RAISED  TO  THE  POWER  OF  TIME  WITH  A 
3 MAX  I  MUM  VALUE  OF  *,F8.1) 

COMPLEX  RZ { 2049 ) ,  FWI2049) 

COMMON  / $ TO R E 1 /  RZ ,  FW 

COMMON  / STORE 2/  FREQ  I  2051  ),  AMPL(?051),  SSI  4098),  S  II 4098) 
C0MM0N/ST0RE3/I P,  YPOS,  YTOT,  TAU,  TINCH,  TPLOT , NAP , TST ( 4 ) , TEND { 4 3 
DIMENSION  SPARE  (4098) 

EQUIVALEE CE  C RZ  C  1 )  , SPARE! 1 )  ) 

READ  (  5,3  000  )  L0G2N,  NAP,  I  PULS,  10GEST,  ZEROLO,  FLOCUT,  HI  CUT, 

1  DELAY,  FNYQ ,  PERIOD 

R  E  AC  (5,999)  TINCH,  TPLOT,  BETA,  (MAX,  ( { TST { I ) , TEND (  I )  )  , I  -  1,4) 


C 

C  1.0G2M 

C 

C  NAP 
C 

C  I  PUL S 

c 

c 

C  LCGEST 
C 

C  ZEROLO 

C  F  LCCUT 
C  HICUT 


2**L0G2N  GIVES  THE  NO.  OF  POINTS  IN  THE  INPUT  AND  OUTPUT 
TIME  SERIES.  ( L  0G2N . LE . 12) 

NO.  OF  TIME  RANGES  FOR  COMPUTING  AUTOPOWER  SPECTRA. 

IF  Is  AP=  0 ,  NO  SPECTRA  ARE  COMIUTED. 

=  1  AND  THE  GRAM-CHARL  IER  WAN  EL  E  T  IS  COMPUTED 

=  2  AND  THE  PULSE  FREQUENCY  RESPONSE  IS  UNITY  FROM  FLf CUT  TO 

HICIT  DECREASING  TO  ZERO  FROf  FLOCUT  TO  ZEROLO 

2*41 OGEST  GIVES  THE  NO.  OF  FREQUENCY  ESTIMATES  FOR  THF 

AUT( POWER  SPECTRA  OF  THE  SYNTHETIC  SEISMOGRAM 

AMP 1  I TUDE  IS  ZERO  FROM  0  TO  \  EROLO  FREQUENCY 

LOW  CUT-OFF  FREQUENCY  FOR  SPICIFYING  SPECTRUM  OF  INC.  PULSE 
H I G(  CUT-OFF  FREQUENCY  FOR  SIECIFYING  SPECTRUM  OF  INC.  PULSE 


- 


ooo  ooo  ooo  o  o  o  o  o  ooo 
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C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


DELAY  TIME  DELAY  TO  POSITION  MAX.  AMPLITUDE  OF  INPUT  PULSE 

FNYQ  THE  NYQUIST  FREQUENCY  FOR  THE  SPECTRA  OF  THE  DATA  TO  BE  CALC. 

PERIOD  THE  APPROXIMATE  PERIOD  OF  THE  INPUT  PULSE  FOR  THE  S.S. 

ACTUALLY  CENTRAL  FREQ.  IS  ABOUT  20%  LESS  THAN  { 1/PERIOD). 

BETA  GAIN  CONTROL  FACTOR  IS  BETA  RAISED  TO  THE  POWER  OF 
TIME  WITH  A  MAXIMUM  VALUE  OF  GMAX. 

GMAX  MAXIMUM  GAIN  VALUE  BY  WHICH  S.S.  IS  MULTIPLIED. 

TINCH  NO.  OF  INCHES  OF  PAPER  PER  SECOND  FOR  PLOTTING  TIME  SERIFS 
T  PLOT  TIME  LENGTH  OF  PLOT  (STARTS  AT  0  5  FOR  SYNTHETIC  SEISMOGRAM 
T ST  START  TIMES  FOR  CALCULATING  AUTOPOWER  SPECTRA  OF  S„  S. 

TEND  END  TIMES  FOR  CALCULATING  AUTOPOWER  SPECTRA  OF  S.  S. 

N  =  2**LCG2N 
M  =  N/2 

MPLUS1  =  M  +  1 
NEST  =  2**LOGEST 
TAU  =  1*/(2.*FNYQ) 

WRITE  (6,1003)  N  ?  M,  FNYQ,  TAU,  PERIOD,  TINCH,  T  PLOT 
WRITE  (  6,1007)  ZEROLO,  FLQCUT  ?  HI  CUT,  BETA,  GMAX 


GENERATE  THE  INPUT  PULSE  AND  TAKE  ITS  FOURIER  TRANSFORM 


CALL  PULSE  { LQG2N, PER  IOD,  IPULS, HICU7, f-LOCUT, DELAY, ZEROLO) 

GENERATE  THE  REFLECTION  COEFFICIENTS  AS  FUNCTIONS  OF  FREQUENCY. 

THIS  REQUIRES  READING  AS  INPUT  VELOCITY  AND  DEPTH  VALUES  PLUS 
Q  AND  Q-DEPTH  VALUES. 

CALL  COEFF  ( L0G2N, FNYQ ) 

MULTIPLY  THE  PULSE  TRANSFORM  AND  COEFFICIENT  VALUES  IN  THE  F  DOMAIN. 


DO  20  1=1 , MPLUS1 
20  FW ( I )  =  F  W (  I )*RZ( I ) 

PUT  FW  IN  SS  IN  A  FORM  TO  USE  FAST  F.T. 

DO  30  1=1, MPLUS 1 
30  SSf  I  )  =  REAL  ( FW ( I ) ) 

DO  40  1=2, M 
IA  =  N  +  2  -  I 

40  SSI  I  A)  =  -  A I  MAG  I FW {  I  )  ) 

PERFORM  THE  INVERSE  FOURIER  TRANSFORM 

CALL  ONE  C  FT  (LOG2N,SS) 

FLN  =  FLCAT(N) 

DO  43  1=1 , N 

43  SP  AR  E (  I  )  =  SSC I ) /FLN 
SI(1)  =  1.0 


. 


n  n  on 
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SS(  1)  =  SS C  n/FLN 
on  50  .1  =  2,  N 

IF  ISI(I-l) .GE.GMAX)  GO  TO  45 


Si  (I  ) 
GO  TO 

=  BETA**{( 1-1 5*TAU) 

50 

45 

SKI) 

=  CM  AX 

50 

ssm 

=  SS(I) *SI ( I ) / FLN 

WRITE 

(6 , 1006) 

TAU 

WRITE 

WRITE 

(6 ,1002 ) 
(6,1001) 

( SI ( I  ) , 1  =  1 ,N  ,2 ) 

WRITE 

(6,1002  ) 

(SS ( I)  , 1  =  1, N) 

PLOT  THE 

SYNTHETIC 

SE ISMOGRAM 

CALL 

PLOT  SS  (4 

, L0G2N , NSAMP ) 

CALCULATE  AUTOPOWER  SPECTRA  OF  SPECIFIED  TIME  INTERVALS 
IF  ( NAP • EQ  *  0 )  GO  TO  100 
DO  53  1=1, N 
53  SSI  I  )  =  SPARE (  i  ) 

LOOP 1  =  LOGEST  +  1 
NES  =  NEST  -5-  1 
DO  55  1=1, NES 

55  FREQ  ( I )  =  FNYQ* U-l ) /NEST 
F  R  E  Q ( N  E  S  T  4  2 )  =  0.0 

FREQ  (  NES  T 1-3  )  =  FLO  AT  I  5* IF  I X ( 1 . / ( 82 . 0*TAU )  +  l) ) 

DO  80  1=1, NAP 

NS T  =  I  F  IX  (  T  ST  (  I  )  /  TAU  )  4-  1 
NEND  =  I  FI  X  {  T END  (  I  )  /TAU  )  4-  I 
NSAMP  =  NEND  -  NST  +  1 
DO  60  I A  =  NST , NEND 
IB  =  I A  -  NST  4-  1 
60  SI ( IB)  =  S  S (  I  A) 

CALL  TAPER  (NSAMP, N) 

CALL  ONE  R  FT  (L0G2N,SI) 

00  7  C  I A = i , N 
70  SPARE  I IA  )  =  SI (  IA) 

CALL  TWO  C  A  C  ( L0G2N , S I  * SP ARE ) 

CALL  PARZEN  (NEST, SI) 

DO  75  I.A--2 ?  NEST 
SI (  I  A)  =  2.0*SI ( IA) 

IB  =  I A  4  NEST- 
75  SI (  IB)  =  0.0 
SI (NES)  =  0.0 

CALL  ONE  R  FT  (LOGPl,SI) 

WRITE  (6,10055  TST(I),  T END ( I  ) 

WRITE  16  1004)  ( ( FREQ ( I  A)  ,  S  I  ( I  A) )  IA=1,NES) 

C 

C  PLOT  THE  AU'OPOWER  SPECTRUM 
CALL  PL0‘  SS  «  5 , LOGP 1  ,  I  ) 

80  CONTINUE 


no  on 


100 
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CONTINUE 

STOP 

END 


SUBROUTINE  PULSE  ( L0G2N, PERIOD, I  PUL S, HI  CUT, FLOCUT ,DELAY , ZEROLO) 

C 

C  THIS  SUBROUTINE  CALCULATES  AN  INPUT  PULSE  FOR  THE  SYNTHETIC  SEISMO- 
C  GRAM  AND  TAKES  ITS  FOURIER  TRANSFORM  UTILIZING  THE  FAST  F.T.  THF 
C  PULSE  IS  A  GRAM-CHARLIER  WAVELET  CALCULATED  FROM  THE  5TH  AND  6 TH 
C  DERIVATIVES  OF  THE  GAUSSIAN  ERROR  FUNCTION*  ALTERNATIVELY,  THE 
C  PULSE  IS  AN  APPROXIMATION  TO  A  DELTA  FUNCTION  AND  IS  SPECIFIED 
C  BY  ITS  AMPLITUDE  SPECTRUM. 

C 

COMPLEX  RZ ( 2  049 ) ,  FWC2049) 

COMMON  / STORE  1/  RZ ,  FW 

COMMON  /  STORE2  /  FREQI2051),  AMPU2051  ),  SS ( 4093  )  ,  SI  (4098) 

COMMON/ ST0RE3/ I P  »  YPOS,  YTGT,  TAU,  TINCH,  T PLOT , NAP , TST { 4 ) , TEND  I  4 } 

1000  FORMAT  { 6 ( 2 X , E 1 4 . 7 ) ) 

1001  FORMAT  {  1H1 ,4X,9HREAL  S I W ) , 8X , 9HI MAG  S ( W) , 8X , 9HAMPL I TUDE , 10X , 5HPHA 
1SE,  1 IX, 9  FFREQUENCY ) 

1002  FORMAT  ('0  THE  AMPLITUDE  SPECTRUM  OF  THF  INCIDENT  PULSE  HAS  BEEN  S 
IP  EC  I F I  ED  *  ) 

1003  FORMAT  Ml  THE  FIRST  300  AMPLITUDE  VALUES  FOR  THE  PULSE  RESULTING 
1  FROM  A  UNITY  SPECTRUM  FOLLOW*//) 

1004  FORMAT  {  10 ( 2X  ,  E  10. 3 )  ) 

1010  FORMAT  (  1H0,53X, *X* , 15X, 1  PULSE!  I  )  5 ,5X,' I f  ) 

1011  FORMAT  (40X,2(5XfE15.7)  » 110) 

N  =  2**L  CG2N 
M  =  N/2 

MPLUS1  =  M  +  1 
FINTL  =  1  ./{ 2.*TAU*M) 

GO  TO  { 5  i  50  )  ?  I  PULS 

A  GRAM  CHARIIER  WAVELET  IS  TO  BE  CALCULATED 
5  TWO  PI  =  6.2831854 
A  =  1 ./SCRT ( TWO  PI ) 

NS  AM  P  =  IFIXl  6MPERI0D/TAU)  +  1 
DO  10  1=1,4098 

io  ssm  =  c.o 

WRITE  (6,1010) 

DO  30  1=)  » NSAMP 

X  =  2.0* I FLOAT! I-1)*TAU  -  3*PER IOP) /PERIOD 

FX  =  ( EX  1  ( - X * X / 2 ) ) *(  ( (  (  (  (-X-M  )*X+i5)*X-10)*X-45)*X+15)*X+15}*A 
SS ( I )  =  ! X 

30  WRITE  (6,1011)  X,  FX,  I 
NSAMP  =  ;*NSAMP 

PLOT  THE  GRJM  CHARLIER  WAVELET 


o  o  o  o  o  o  o 
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CALL  PLOT  SS  I  1  ,  LOG2N , NS AMP ) 

CALL  ONE  R  FT  (LOG2N,SS) 

SS  NOW  CONTAINS  THE  COMPACT  STORED  FOURIER  TRANSFORM  OF  THE 
INPUT  PULSE 

CALL  RECCVR  { L OG2N  ,  S  S  ,  S I ) 

THE  REAL  PART  OF  THE  F ,T .  IS  NOW  IN  SS,  THE  I  MAG.  PART  IN  SI 

WRITE  ( 6  s 1 00  1 ) 

DO  40  I=1,MPLUS1 
FREQ ( I )  =  i  I-1)*FINTL 
F  W (  I  )  =  CMP  L  X ( SSII),SI(I)) 

AMPL ( I )  =  CABS ( F  W (  I )  ) 

IF  ( (SI ( I) .EQ.O.O) .AND. I SSII) .EQ.O.O) )  SS  f  I )  -  1.0 
PHASE  -  AT  AN2 ( S I (I) , SSI  I) ) 

IF  { I . GT  .200)  GO  TO  40 

WRITE  (6,1000)  FW(I) ,  AMPL(I),  PHASE,  FREQ { I ) 

40  CONTINUE 

C  PLOT  THE  AMPLITUDE  SPECTRUM  OF  THE  GRAM  CHARLIER  WAVELET 
CALL  PLOT  SS  I  2 , L0G2N ,NSAMP ) 

RETURN 

C 

C  A  DELTA  FUNCTION  INPUT  PULSE  IS  TO  BE  CALCULATED.  THE  AMPLITUDE 
C  CHARACTERISTICS  ARE  FIRST  COMPUTED. 

50  LZ  =  I F  I  X ( ZERO LG /F INTL) 

LY  =  LZ  +  1 

IF  ILZ.EC.O)  GO  TO  52 
DO  51  1=1, LZ 
FW II)  =  (0. 0,0.0) 

AMPL (  I  )  -  0.0 

51  SS ( I )  =  0.0 

52  LA  =  IFIXIELOCUT/F INTL )  +  1 
LB  =  LA  +  1 

LC  =  IF  I  X(  H I  CUT /F  INTL  )  +  LA  -  LY  i-  1 

LO  =  LC  -  LA  +  LY  ~  I 

LE  =  LC  +  1 

PI  =  3.14159265 

ARG  =  0. 5*PI/I LA-LY) 

DO  55  I  =  .  Y  ,  LA 
IB  =  I  -  LY 
IA  =  LC  -  IB 
SINE  =  S  INI ARG* IB) 

FWC I  A)  =  CMPLX( S INF, 0 .0) 

AMPL { I )  =  SINE 

SSI  I)  =  SINE 

FWII)  =  :MP LX  I  SINE, 0.0) 

AMPL ( I  A )  =  SINE 


■ 


■) 


. 


no  on  on 


55  S  S ( I  A }  =  SINE 
DO  60  I  =  LB  t  L  D 
F  W  {  I  )  =  11.0,0.0) 

AMPL  ( I )  =  1.0 
60  SS ( I )  =  1 .0 

DO  65  I  =  LE , MPLUS 1 
FW ( I )  =  (0,0,0,01 
A  M  P  L  (  I  )  =  0,0 
65  SS ( I )  =  C . 0 

DO  68  1  =  NPL  US  I ?  N 
68  SS (  I  )  =  0.0 
C 

C  AN  INVERSE  F.T.  PRODUCES  A  PULSE  IN  THE  TIME  DOMAIN 
CALL  ONE  C  FT  (L0G2N, SS) 

FLN  =  FLOAT  IN) 

DO  70  1=1, N 
70  SS {  I)  =  SSI  I ) /FLN 

A  DELAY  IN  TIME  OF  THE  PEAK  AMPLITUDE  OF  THE  PULSE  IS  MADE 
IF  (DELAY. EQ. 0.0)  GO  TO  85 
LDEL  =  I F I X { DEL  AY/TAU )  +  1 
LDEL1  =  LDEL  +  1 
DO  75  1=1, LDEL 
IA  =  LDEL  -1+1 
75  SI(  I)  =  SS  C  I  A) 

DO  80  I  =  L  DEL  1 , N 
IA  =  I  -  LDEL 1  +  2 
80  SICI)  =  SS(IA) 

DO  83  1=1, N 
83  SSI  I )  =  SI {  I  ) 

85  WRITE  (6,1003) 

WRITE  (6,1004)  { SS ( I  ) ,  1  =  1  ,  500) 

NS AMP  =  L  . 0 / { T I N C H *T  A U ) 

PLOT  THE  INIUT  PULSE 

CALL  PL 01  SS  ( L  ,  L0G2N  ,  NSAMP ) 

THE  COMPLEX  SPECTRUM  OF  THE  PULSE  IS  COMPUTED. 

DO  100  1=  1 , MPLUS1 
100  FREQ ( I  )  -  (  I-I ) *FI NTL 

WRITE  (6-1002) 

IF  (DEL A'  .EQ.O.O)  GO  TO  120 
CALL  ONE  R  FT  (L0G2N,SS) 

CALL  RFCfVR  (L0G2N , SS , SI ) 

WRITE  (6,1001) 

DO  110  1 :  1 , MPLUS 1 
FW ( I )  =  (  MPLX(SS( I >  T  S I (  I  )  ) 

AMPL (  I  )  CABS  (  FW (  I  )  ) 

IF  (  (SI  (  ) .EQ.O.O) .AND. ( SSI  I  ) .EQ.<  .0)  )  SS(I)  =  1.0 
PHASE  =  / T  AN2 I S  I  ( I ) , S S (  I )  ) 


n  n  n  n 
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IF  (  I  .GT  .200 )  GG  TO  110 

WRITF:  (6,1000)  FW(I),  AMPL(I),  PHASE,  FREQ{  I  } 
110  CONTINUE 


C 

C  PLOT  THE  AMPLITUDE  SPECTRUM  OF  THE  INPUT  PULSE 
120  CALL  PLOT  SS  12 » LOG2N ,NSA MP ) 

RETURN 

END 


SUBROUTINE  COEFF  { L0G2N , FNY Q ) 

THIS  SUBROUTINE  GENERATES  THE  REFLECTION  COEFFICIENTS  OF  THE 
LAYERED  MEDIA,  WITH  ATTENUATION,  AS  A  FUNCTION  OF  FREQUENCY. 


COMPLEX  RZ I  2049 } 5  FW 12049) 

COMMON  / STORE  1 /  RZ ,  FW 

COMMON  /  ST0RE2/  FREQ(2051),  AMPLI2051),  SSI4098),  S  It  4098) 

COMMON/ S TOR E3/ I P ,  YPOS,  Y  TOT ,  TAU,  TINCH,  TPLOT , NAP , TST ( 4 ) , TEND ( 4 ) 
REAL  VR ( 600 ) ,  Z(600),  Q(20),  ZQ(20) 

COMPLEX  E 1 600 ) ,  VEL(  600),  S ( 600 ) ,  R,  D,  F,  E,  RZERO 

1000  FORMAT  (215) 

1001  FORMAT  (  lX,E14.7,iX,E14.7,llX,E14.7,lX,E14.7) 

1002  FORMAT  { £ F  1 0  •  0 ) 


1008  FORMAT 

1009  FORMAT 

1010  FORMAT 

1011  FORMAT 

1012  FORMAT 
1 » PHASE 

1013  FORMAT 

1014  FORMAT 
N  =  2**L(G2N 
M  =  N/2 
MPLUS1  =  M  + 
READ  (5,  -COO) 
NZ1  =  NZ  +  1 
READ  (5f:001) 
READ  ( 5,  002 ) 


(*0  THE  COMPLEX 


SLOPE 


VALUES  FOLLOW 


(4 


w  a  i  line  /i 

V  Ml.  U  O  /  l_. 


I N 


c  ) 


{ 4( 4X,E12 
(  *  1  DEPTH 
{ 2 ( 4X , E 1 2 
(  1H1 , 24 X , 

OF  R  8  , 4X , 

( 21X , 5E16 .6 ) 

{ ' 0  Q-DEPTH  VALUES 


4,2X,E12.4) ) 

VALUES  AND  COMPLEX 

4 


VELOCITIES 
4 , 6  X  ,  E  1 2 . 4 , 2  X  ,  E  l  2 . 4 )  ) 

'FREQUENCY8  ,8X, ‘REFLECT  ION  COEFF I  Cl  ENT* 
‘AMPLITUDE  OF  R 1 ) 


{ 2  SETS/LINE ) 5  ) 


8X 


AND  Q  V/LUES  (4  SETS/LINE)') 


1 


C 

c 

c 

c 

c 

c 

c 

c 


NZ 

NZQ 

VR  (  J) 
Z  ( J  ) 
Q(  J) 
ZQ  ( J) 


(  5, 

THE 

THE 

THE 

AND 

THE 

AND 


NZ,  NZQ 

( VR ( J ) , Z ( J ) 
{ Q ( J ) , ZQ  (  J  ) 


J-2 , N  Z  j  ) 
J= 1 »  NZO ) 


AND 


DEPT  I 
V  AH 


NO.  OF  VELOCITY 
NO.  OF  Q  AND  Q-DEPTH 
VELOCITY  VALUES 

THE  CORRESPONDING  DEPTH  VALUES 
Q  (INVERSE  ATTENUATION)  VALUES 
THE  CORRESPONDING  Q-DEPTH  VALUES 


VALUES  TO  BE  READ 
ES  TO  BE  READ  IN 


IN 


VR(  1)  =  ''R  (2) 
Z  (  1  )  =  0 


i  T 


*. 
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IA  =  1 

DO  iO  1=1, NZ1 

IF  I  Z { 1  )  cLE  *  ZQ{  IA)  }  GO  TO  5 
I A  =  I A  +  1 

5  VI  =  0.5*VR(  I  )  /  QUA) 

10  VE  L ( I )  =  CMPLX( VR(  I)  , VI  ) 

S ( 1 )  =  CMPLX( 0.001 ,0.05 
S  ( N  Z 1 )  =  S  {  1  ) 

WRITE  16,1010) 

WRITE  16,1011)  UZU)  ,  VEL  (  I  )  )  » I =1 , NZ I ) 

WRITE  (6,1014) 

WRITE  (6,1009)  { (ZQ{ I  )  ,Q(  I)  )  ,I  =  1,NZQ) 

DO  20  1=2, NZ 

SCI)  =  (VELCI  +  1)  -  VEL  (  I  )  )  /  (Z(I+1)  ~  ZU)) 

SR  =  REAL ( S ( I )  ) 

IF  ( ABS ( SR)  .GE. 0.001 )  GO  TO  20 
SI  I  =  A I  M  A  G  (  S  (  I  )  ) 

SCI)  =  CMPLX ( .001 , SI  I  ) 

20  CONTINUE 

WRITE  (6,1008) 

WRITE  (6,1009)  ( S< I ) , 1=1, NZ1) 

WRITE  (6,1012) 

DO  100  1=1, M PLUS! 

W  =  6. 283185 4*FREQ( I  ) 

IF  ( I  . EQ  -  1 )  GO  TO  60 
DO  30  I A  =  1 , N Z  1 

30  B ( I A )  =  CSQRTCl.  -  4 . * ( W/ S (  I  A ) ) **2 ) 

R  =  (0  *0  ,0.0  ) 

DO  50  I  A  =  1 , N Z 
IB=NZ-IA+  1 
IC  =  IB  +  1 
D  =  SC IB)*B(  IB) 

E  =  S(IC)  -  S  (  I  B ) 

F  =  S(IC)  *B(  IC  )*(  l.-R)  /(l.  +  R) 

IF  ( (REALi VEL(  IC ) ) .EQ.REALl VEL (  IB)  ) ) . AND . ( A IM AG ( VEL ( IC)  ) .EQ. 
1  A  IMAG  ( VE,_  (  13)  )  )  )  GO  TO  40 

R  =  ((D+E-F)  /  (D-E+F))  *  CEX P { -B ( I B ) *CLOG ( VEL ( I C ) / VEL ( I B ) ) ) 
GO  TO  50 

40  R  =  (D+E-F)  /  (D-E+F) 

50  CONTINUE 
RZERO  =  -I 
GO  TO  70 

60  RZERO  =  ( VEL  ( 1  )  -  VEUNZD)  /  (VEL(l)  +  VEL(NZ.l)) 

70  RZ(  I  )  =  RZERO 

AM  PL ( I )  =  CABS (RZERO) 

PHASE  =  A  f AN2  ( A IM AG ( RZERO)  , RE AL ( R ZERO ) ) 

FR  =  W/6. 2831354 

WRITE  16,1013)  FR,  RZERO,  PHASE,  AMPL(I) 

100  CONTINUE 


C 


noon 
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C  PLOT  THE  TRANSFER  FUNCTION  OF  THE  LAYERED  MEDIUM 
CALL  PLOT  SS  ( 3  ,  LOG2N  ,  NSA MP ) 

RETURN 

END 


SUBROUTINE  PLOT  SS  ( I  TYPE ,LOG2N , NSAMP ) 

THIS  SUBROUTINE  PLOTS  THE  VARIOUS  QUANTITIES  CALCULATED  ELSEWHERE, 
INCLUDING  THE  FINAL  SYNTHETIC  SEISMOGRAM. 

1000  FORMAT  (c0  THE  TOTAL  NUMBER  OF  CALCOMP  PLOTS  FOR  THE  PLOTTER  REQUE 
1ST  SLIP  IS  1 , 15) 

DIMENSION  BUFF { 2048 ) 

COMPLEX  RZ ( 2049 ) ,  FWC2049) 

COMMON  / STORE  1 /  RZ ,  FW 

COMMON  /  ST0RE2/  FREQ(2051  ),  AMPU  2051  ),  SSI4098),  SK4098) 

COMMON/ STORE 3  /  I  P  ,  YPOS,  Y  TOT ,  TAU,  TINCH,  T PLOT , NAP , TST ( 4) , TEND { 4 ) 
N  =  2**LCG2N 
M  =  N/2 

MPLUS1  =  M  +  1 

MPLUS2  =  M  +  2 

M PL US 3  =  M  +  3 

GO  TO  (20,  40,  60,  80,  120),  ITYPE 
20  CALL  PLOTS  (  BUFF ( 1 ) , 8 192 ) 

IP  =  1 
YTOT  =0.0 

CALL  PLOT  (2.0, 1.0, -3) 

IP  =  IP  +  1 
NS  1  =  NSAMP  +  1 
NS2  =  NSAMP  +  2 
DO  25  1=1, NSAMP 
A  MPL (  I )  =  SS  (  I  ) 

25  FREQ  (  I  )  =  FLOATU-1 )  *TAU 
FREQ(NSl)  =  0.0 
FREQ ( NS2 )  =  1. /TINCH 
CALL  SCALE  ( AMPL ,4 .0, NSAMP, 1 , 10 . ) 

CALL  AXIS  (0.0, 0.0, 'TIME  ( S EC ONDS ) !  , - 1 4 1 4 . 0 , 0 . 0 , FRE Q ( NS  1) , 
1FREQ(NS2 ) , 10. ) 

CALL  AXIS  ( 0 .0,0.0 ,» PULSE  AMPLITUDE  1 , 15 ,4 . 0, 90 . 0 , AMPL ( NS1 ) , 

1 AMPL ( NS2 ) i 10. ) 

CALL  LIME  (FREQ, AMPL, NSAMP, 1,0,3) 

CALL  SYMBCL  ( 0 . G  ,4 . 5 , C . 25  , '  I NPUT  P l L S E *  ,  0 . 0 , 1 1 ) 

YPOS  =  8.C 

YTOT  =  YTOT  +  YPOS 

RETURN 

40  CALL  PLOT  ( 0 . 0 , Y POS , -3 ) 

IP  =  IP  +  1 

FREQ  ( MPLLS2 )  =  0.0 
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FREQ  ( MPLUS3 )  =  FLOAT (5  *  IF  IX (1 . / ( 82  * 0*TAU) +1 ) ) 

CALL  SCALE  ( AMPL , 6 . 0 , MPLUS1 , 1 , 10 . ) 

CALL  AXIS  (0.0 ,0.0, 5  FREQUENCY ( HZ . )* , -  15 , 8 . 0 , 0 . 0 , FRF Q( MPLUS2 ) , 

1  FREQ(MPLUS3) ,1C.) 

CALL  AXIS  ( 0 .0 ,0 .0  ,* PULSE  TRANSFORM  AMPLITUDE », 25 ,6 . 0, 90 . 0, 

1  AMPL ( MPLUS2 ) ,AMPL(MPLUS3 ) ,10.) 

CALL  LINE  (FREQ,  AMPL,  MPLUS1 ,  1,  0,  3) 

CALL  SYMBOL  (0.0,6.5,0.25, ’AMPLITUDE  SPECTRUM  OF  INPUT  PULSE*, 

1  0.0,33) 

YPOS  =  1C.0 

YTOT  =  Y  TOT  +  YPOS 

RETURN 

60  CALL  PLOT  ( 0 . 0 , YPOS ,-3 ) 

IP  =  IP  +  1 

CALL  SCALE  ( AMPL ,6 . 0, MPLUS1 , l , 10. ) 

CALL  AXIS  (0.0, 0,0,*  FREQUENCY ( HZ  .  ){  ,- 15 , 8 . 0, 0 . 0, FREQ( MPLUS2 ) ? 

1  FREQ ( M  PLUS3 )  ,10.) 

CALL  AXIS  { 0.0 ,0.0 ,* AMPLITUDE  * , 9 ,6 . 0, 90. 0, AMPL (MPLUS2 ) , 

I  AMPL(MPLUSB) , 10. ) 

CALL  LINE  (FREQ,  AMPL,  MPLUS1,  1,  0,  3) 

CALL  SYMBOL  (0.0,6.5,0.2  5, *  TRANSFER  FUNCTION  OF  LAYERED  MEDIA’, 

1  0.0,34) 

YPOS  =  -  YTOT 
YTOT  =  0. 

RETURN 

80  CALL  PLOT  { 1 1 . 0 , YPOS, -3 ) 

IP  =  IP  +  1 

N  =  IF  IX  (TPLOT/TAU  )  -8-10 
IF  (N.GT.2**L0G2N)  N  =  2**LOG2N 
DO  83  1=1, N 
83  SI  (  I )  =  SS(  I  ) 

CALL  SCALE  ( SS  ,  4 .0  ,N, 1 , 10  . ) 

SS(N-K)  =  -  SS  (  N  +  l ) 

NPLUS2  =  N  +  2 

SS ( NPLUS  Z )  =  -  S  S { NPLUS2 ) 

DO  85  1=1, N 

85  SS(  I)  =  S S (  I)/SS(NPLUS2) 

T I  NCR  =  1 .0 /T INCH 

YD  I  ST  =  T PLOT *T INCH 

NLINES  =  IF  I  X ( YDI ST/25.05  )  +  1 

XP  =  0.0 

T2  =  0.0 

T 1  =  0.0 

N2  =  I 

T  =  TINCH  *  TAU 
DO  110  I -1, NLINES 
IF  ( YDIS1 .LE.25.0)  GO  TO  90 
YD  =  25.0 

YD  I  ST  =  YDIST  -  25.0 
GO  TO  95 


90  YD  =  YD I  ST 
95  T 1  =  T 1  +  12 

12  =  T2  +  YD/TINCH 

CALL  AXIS  IXP,0.»'TIME  ( SECONDS  )  5  » 14»YD,90.0,T1,TINCR» 10. ) 
CALL  AX  I  S(  XP+.  5 ,0*  0,  *  OUTPUT  AMPLITUDE'  ,-16,4.0,0.0,SS(N*l), 

1  SS  C.NPLUS2  )  ,10.  ) 

CALL  PLOT  ( XP+2.5,  0. *  -3) 

IP  =  IP  +  1 
XP  ■=  3.0 
N 1  =  N2 

N2  -  NX  +  I F I  X C I T2-T1 J/TAU) 

DO  100  I  A=N 1 f N2 
TPOS  =  T*FLOATl IA-N1 } 

100  CALL  PLOT  ISSCIA),  TPOS,  2) 

110  CONTINUE 

CALL  PLOT  ( XP , 0 . 0 , —3 ) 

IP  =  IP  4  1 

CALL  SYMBOL  { 0 • 0 , 0 . 0 , 0 . 3 0 , * SYNT HET I C  SEISMOGRAM® ,90.0,20) 

CALL  PLOT  (4  .0,0. 0,-3) 

IP  =  IP  +  1 
DO  115  I  - 1  ,  N 
115  SS(I)  =  SI(I) 

IF  (NAP.NE.O)  RETURN 
CALL  PLOT  (0.0,0.0,999) 

ViRITE  (6  i.IOOO)  IP 
RETURN 

120  IE  (NSAMP.NE.l)  GO  TO  130 

CALL  AXIS  10.0,0.0,' FREQUENCY  (  HZ  .  )  *  , -15 , 8. 0, 0. 0 , FREQ ( MPLUS2  ) 
1  FREQ ( MPLUS3 ) , 10. ) 

SMIN  =  SKI  ) 

DO  12  2  I -2 , MPLU  S 1 

IF  {  SMIN.GT.Sn  I  )  )  SMIN=SI(I) 

122  CONTINUE 

DO  125  L  It MPLU SI 
IF  (SIC  I  .EQ.O.O)  SHI)  =  SMIN 
125  Sill)  =  /  LOG  IOC SI(I)5 
YPO  =  OK 
XPO  =  0.  C. 

GO  TO  140 

130  IFCNSAMP  NE.2)  GO  TO  133 
YPO  =  10  0 

XPO  =  OJi 

YPOS  =  0  5/ (TAU*MPLUS1 ) 

N  1  =  I  F  I  ;  (  5-.0/YPOS)  +  1 
N2  =  I  F  I  :•  (25.0/YP0S) 

DO  131  I-Nl , N2 
IA  =  I - N  l  +  l 

131  FREQ ( IA)  =  FREQ l  I  ) 

MPLUS1  =  N2  -  N1  +  1 
MPLUS2  =  MPLUS1  +  1 
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MPLUS3  =  MPLUSl  +  2 
FREQ  (  MPL  LS2  )  -=5.0 
FREQ ( MPL IS3 )  =  2.5 
GO  TC  135 

133  YPO  =  2c C 

XPO  =  -0.7*(NSAMP-2) 

N1  -  IFI X( 5.0/YP0S)  +  1 
N2  =  IFIX(25.0/YPOS) 

MPLUSl  =  N2  -  N1  +  1 
NPLUS2  =  MPLUSl  4  1 
MPLUS3  =  MPLUSl  4-  2 
135  Mil  =  N1  +  1 
SHIN  -  SI  (Nl) 

DO  134  I - N 1 1 ,  N  2 

IF  (SMIN. GT. SI(I)  )  SMIN  =  SHI) 

134  CONTINUE 

DO  137  I =N1 , N2 

IF  (  SI  (  I  )  .  EQ . 0 «  C )  Sill)  =  SMIN 
IA  =  I-Nl+1 

137  SI  (  I A  3  =  AL0G10  (Sim) 

140  CALL  PLOT  (0«0,YPG,~3) 

IP  =  IP  +  1 

IF  (NSAMP.EQ.2)  CALL  AXIS  I  0 .0 , 0 . 0 ,» FREQUENCY  ( HZ . ) 1 ,-15 , 8. 0, 0. 0 , 

1  FREQ l MPL US 2)  ,  FREQ ( M PLUS 3 3,10.) 

CALL  SCALE  ( SI , 8 .0 ,MPLUS 1 , 1 , 10. ) 

CALL  AXIS  {XPOtO.O, * AUTOPGWER  SPECTRAL  DENSITY  l LOG  10  SCALE)*, 41, 
18. 0,90.0  ,$I ( MPLUS2 )  , S I ( MPLUS3) , 10.0) 

CALL  LINE  ( FREQ, SI , MPLUSl, 1,10, NSAMP) 

YY  =  I  SI (MPLUSl 3  -  SI ( MPLUS2 )  )  /  SIIMPLUS3) 

ST  =  TST(NSAMP) 

SE  =  TEND l NS AMP ) 

CALL  NUMBER  (8.5,YY,.15,ST,0.0,2) 

CALL  WHEI E  l XX, YZ, FCTR) 

CALL  S.YMLOL  {XX,YY,.15,'  TO  *,0.0.4) 

CALL  WHE  E  ( XX, YZ, FCTR) 

XX  =  XX  0.1 

CALL  NUM'IER  { XX  ,  YY  ,  .  1  5  ,  S E  , 0  . 0 , 2  ) 

CALL  WHERE  (XX, YZ, FCTR) 

CALL  SYMTOL  (  XX  »  YY  ,  .  1.  5  ,  *  SEC.»,0.nf5) 

IF  (NSAM 5.NE. NAP)  RETURN 
IF  l  NAP.  vE.  1  )  GO  TO  150 
YPO  =  0.  ) 

GO  TO  163 

150  YPO  =  —2 .0*1  NAP-2)  -  10.0 
160  CALL  PLOT  ( 1 1  .0  , YPO , -3 ) 

IP  =  IP  M 

CALL  PLOT  (0.0,0.0,999) 

WRITE  (6,1000)  IP 

RETURN 

END 


’ 

* 
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SUBROUTINE  TAPER  (NSAMP,N) 

THIS  SUBROUTINE  TAPERS  THE  FIRST  AND  LAST  10?  OF  THE  DATA  FOR  WHICH 
THE  AUT GROWER  SPECTRA  IS  TO  BE  COMPUTED, 

COMMON  /STORE 2/  FREQ (2051),  AMPLI2051  ),  SSI  4098),  SI  (4098) 

COMMON/STORE3/IP,  YPOS,  Y  TOT,  TAU,  T I NCH  f  TPLOT * NAP , TST ( 4 ) , TEND { 4 ) 

NS  1  =  NS  AMP  4-  1 

NS  10  =  N  SAM  P/ 10 

PI  =  3,  I  A 1592 65 4 

ARG  =  5-C*PI/NSAMP 

DO  10  JA  =  1 ?  NS  10 

JAA  =  NS  AMP  -  JA  f  1 

S  -  SI N  f ARG*{ JA-1)  ) 

SKJA)  =  S*SI(JA) 

10  SI (JAA)  =  S*SII JAA) 

DO  20  JA  =  NS 1 f  N 
20  SKJA)  =  0.0 
RETURN 
END 


SUBROUTINE  PAR  ZEN  (  M,  C  ) 


MODIFIES  THE  COVARIANCE  FUNCTION  C(T)  BY  THE  PARZEN  LAG  WINDOW. 
THE  RESULTS  ARE  STORED  IN  C(T),  T= 1 , 2 , 3 , . . . , M+ 1 . 

REAL  C  (10),  D 

INTEGER  J’t  M  UPON  2,  TWO  M 

M  UPON  2  =  M/2 
DO  10  J  = 1 , M  UPON  2 

D  =  1-0  -  6.0*{  FLOAT ( J ) /FLOAT  M)  )**2 
.  +  6 .  C  *  l  FLGATU5  /FLOAT  M)  )**3 

C{ J  +  l  )  =  D  *  C( JFl  ) 

10  CONTINUE 

M  UPON  2  =  M  UPON  2  4-1 
DO  20  J=1  UPON  2,M 

D  =  2.04(1.0  -  FLOAT!  J  ) /FLOAT  (II)  )**3 
C { J  + 1  l  =  D  *  Cl J  +  l) 

CONTINUE 
R  E  T  U  1  N 
END 


20 


. 
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SUBROUTINE  TWO  C  A  C  l L0G2Nf  X  *  Y } 

C  TWO  CYCLIC  AUTO  COVARIANCES 

C 

INTEGER  L0G2N 
REAL  X  (  1.0)  ,  Y  (10) 

C 

INTEGER  J  f  K  t N  »  N  UPON  2 
REAL  P 1 t  P 2 
C 

N=2**LGG2N 

X(1)  =  X(  1  )**2 
y { i  j  —  v ( i j#«2 
N  UPON  2=N/2 
DO  100  J  =  2  t  N  UPON  2 
K=N+2- J 

P1  =  X ( J )**2+X ( K ) **2 
P2=Y(J)**2+Y(K)**2 
X ( J  )  =  P 1 
X  (  K )  =  P 1 
Y ( J  )  =P2 
Y  <  K  )  =  P  2 
100  CONTINUE 
K=N/  2-t- 1 
X  (  K  )  =  X  (  K  )**2 
Y(K)=Y(K}**2 

CALL  MR  ID  FT  IL0G2N,X,Y) 

RETURN 

END 


SUBROUT  I b E  ONE  R  FT  (L0G2N,X) 

C*************?***************************************************** 

C  ONE  REAL  FOURIER  TRANSFORM 

C 

INTEGER  LGG2N 
REAL  X  (10) 

C 

INTEGER  x  tJNfK,KN,N»N  OVER  2 
REAL  ARG,CtPI ,S,Tf XI,XRfYI,YR 
C 

PI -3*141! 92654 
N=2**( L0G2N-1 ) 

CALL  R  S LB  B  0  (L0G2N?X) 

CALL  R  SUB  B  0  { L0G2N- 1 , X ( 1  )  ) 

CALL  R  SIB  B  0  {  L0G2  N- 1  ,  X  {  N  4-1  )  ) 

CALL  MR  ID  FT  ( L0G2N- 1 , X ( 1 ) , X { N+ 15) 

N  OVER  2--N/2  +  1 


r>  no 
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00  ICO  J  =  2,N  OVER  2 
K=N  +  2- J 
J  N  =  J  +  N 
KN=K+N 

XR=( XI J)  +X( K ) )*.5 

XI=(X( JN)-X(KN) )*.5 

YR=(X{ JNHX(KN) )*.5 

YI=IX(K)-XCJ } )*.5 

ARG  =  P I  AFLOAT { J- 1 ) /FLO AT  I N 3 

C=C0S ( ARG ) 

S=S  INI ARG) 

T=YR*C+Y I*S 
Y I=Y I *C~ YR*S 
YR=T 

X{ J ) =XR+ YR 
X I K }  =  XR-YR 

X  (  KN)  =  X  I  +YI 
XI JN)=YI-XI 

100  CONTINUE 

XR- X { 1 J+XIN+l ) 

YR=X( 1 ) -XI N+l ) 

XI  1  )-=XR 
XI N  +  l  )=YR 
RETURN 
END 


SUBROUTINE  R  SUB  B  0  ( L0G2N t X ) 

REVERSE  SUBSCRIPT  BIT  ORDER 

INTEGER  L0G2N 
REAL  X  (10) 

C 

INTEGER  .J 
REAL  T 
C 

INTEGER  /fB,CtD,E,F,G,H,IfJ»KtLfM,N,BS,CS,DS,ES,FSfGS,HS,IStJS,KSt 
.LS>MSf NS,AL» BL,CLf DLf ELtFL,GLfHL, IL, JLtKLtLL,ML»NLt  S  (14),  U  (14) 
EQUI VALET CE  l BS , S ( 2 ) ) , ( CS » S l 3 ) ) , ( [ S , S ( 4 ) ) , ( ES , S ( 5 ) ) , I FS , S ( M ) , 
•CGS,S(7))t(HS»S(8))t(IS»S(9))»{JSiSI10))j(KS»S(ll))?(LS»S(2))f 
.  (MS*  S(  13  1  )  ,  (NS  tS(  14)  )  ,  (AL,U(  1 )  >  t  IE  L,U(2  )  I  *  (CL  f  UI  3))  »  (DL,U(‘  - )  )  t 
. (EL,U(5  )  1 , { FL,U (6)  ) T  (  GL,U (7) ) * ( HL , U{ 8  )  )  ,  (  IL,U( 9) ) , I JL»U( 10  ) , 

•  IKL,U( 11 )) , I LL,U( 12)  )  , ( ML , U ( 13 ) ) , I NL,U( 14) ) 

C 

NS=2**(L( G2N-1  ) 

NL=2*NS 
DO  100  K~  2 » 13 
J=  1  5-K 


d  n  n 
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U( J1=S( J+l) 

S( J}=1 

IF ( S ( J  +  1 ) *  GT • 1 )  S(J)=S(J+l)/2 
100  CONTINUE 


AL  = 

8S 

JJ  = 

0 

DO 

200 

A  =  1 ,  A  L 

DO 

200 

B=A , B  L , B  S 

DO 

200 

C  =  8 , CL , C  S 

DO 

200 

D=C , DL , D  S 

DO 

200 

E  =  D , EL , E  S 

DO 

200 

F  =  E  »  F  L  ?  F  S 

DO 

200 

G  ~  F  ,  G  L  ,  G  S 

DO 

200 

H=G,HL,HS 

DO 

200 

I=H,IL»IS 

DO 

200 

J=I,JL,JS 

DO 

200 

K=J,KL,KS 

DO 

200 

L=K, LL,LS 

DO 

200 

M=L , ML , MS 

DO 

200 

N-M, NL,NS 

JJ  = 

JJ  +  1 

IF 

(JJ. 

LE.N)  GO 

T=X ( J J  ) 

XI J J)=X( N) 
X ( N  )  =  T 

200  CONTINUE 
RETURN 
END 


SUBROUTINE  MR  ID  FT  (L0G2N,X,Y) 

MIXED  RADIX  ONE  DIMENSIONAL  FCURIE1  TRANSFORM 

INTEGER  L3G2N 
REAL  X  I  13) ,  Y  (10) 

C 

INTEGER  JJ,JO, J1,J2, J3,N,M4 

REAL  ARG  ,C1 ,C2»C3» I0,Ii,I2,I3,R0,Rl,R2,R3,Sl,S2,S3fT 
C 

INTEGER  * ,B, C,D,E,F,G,H, I , J,K,L,M, BS,CS,DS, ES, FS,GS,HS, IS,. S,KS, 
•  LStMSt  ALf3L,CL,DL,EL,FL,GLtHL,ILt  JL,KL,LL,ML,  S  (13),  IJ  (15) 
EQUIVALENCE  CBS,S(2)),(CS,S{3)),(CS,S(4)},(ES,S(5)},(FS,S(M), 

.  (  G  S  »  S  (  7  )  ) , I HS , S ( 8 1  )  ,  (  IS?S(9)),(JS,S(10)  )  , ( K  S , S (  1 1 )  )  ,  ( LS , S ( '  2)  ), 

. (MS,  St  13  )  )  , ( AL,U(  1  )  ) , (BL, U( 2)  ) , (Cl ,U 13)  ) , ( DL , U(4) )  ,  (EL ,U (5  J  ) , 
•(FL,U(6)  )  , ( G  L , U ( 7 )  )»(HL,U{8)),{IL,U(9)),(JL,U(10)),(KL,U(li  )), 

. (LL,U(  12) >  , (ML, U( 13)  ) 

C 

N=2**LOG2N 


. 


100 

200 

300 

400 

500 
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IF  { LQG2N®  L£  *  1 )  GO  TO  500 
CO  400  K=2,LOG2N,2 
M=2**{  LGG2N-K) 

M4=4*M 

DO  300  J=1,M 

ARG=6. 2831 853 07*FLGAT(J-1  ) /FLOAT  ( N4  ) 
C1=C0S( ARG) 

S1  =  S IN { ARG  3 

C2=C1*C1-S1*SI 

S2=C1*S1*C1*S1 

C3=C2*C1-S2*S1 

S3=C2*S1 +S2*C 1 

DO  200  I=M4f Nf M4 

J0-I+J-M4 

J1=J0+M 

J2=J1+M 

J3= J2+M 

R0=  X { J  0  3 +  X( J  2 ) 

R1=X(J0)-X(J2) 

10  =  Y< JO) +Y{ J2) 

11  =  YC J03-Y( J  2 ) 

R2  -  X i J 1 ) +X{ J3) 

R3=XI Jl)-X( J3) 

I2=Y(  Jl)  +YU3) 

I  3  =  Y  {  J  1  )  -  Y  {  J  3  ) 

XI JC)=R0+R2 
Y(JO>=IO+I2 

IF  (  ARG.  EQ.O.O)  GO  TO  3.00 
XI J2)=(R1  +  I3 )*Cl  +  (  I  1-R3) * S 1 

Y  {  J  2  )  -  (  I  1-R3)*C1-(R1+I3)*S1 
X ( J 1 )  =  ( R 0-R2 ) *C2+ ( I  0- 1 2 ) *S2 
Y ( J 1 }  =  U  C- I  2 ) *C 2 - { RO-R  2 ) *  S2 
X(J3)=(R'-X3)*C3+(I1+R3)*S3 

Y  (  J 3 )  =  1 1  4-R3)*C3-t  R1-I3)  *S3 
GO  TO  20 G 

CGNT INUE 
X  {  J  2  )  =  R  1-  13 

Y  {  J  2  )=  1 1-  R3 
X( J1)=R0- R2 

Y  {  J  1 3=  I  0-  12 
X{  J3)=R1  -13 
Y { J3 )  =  I  1  R3 
CONI  INUE 
CONT INUE 
CONTINUE 
CONTINUE 

IF  ( L0G2  oEQ  eLGG2N/2$2 )  GO  TO  700 
DO  600  I  1,N?2 

ro=x { i )+;:( i  +  i) 

Ri=x  ( I i  +  i ) 


o  o  o  o 
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I  0 -V  i  I  KYI  1  +  1) 

I1=Y ( I )-Y( 1+1) 

X( I )=R0 

Y  (  I )  =  I  0 
X{  I  +  1 ) —  R  1 

Y  (  I  +  1 )  =  I  1 
600  CONTINUE 
700  CONTINUE 

MS-N /2 
ML  =  N 

DO  800  K=2,12 
J= 14-K 
S( J5=l 
U( J ) =S { J+l) 

IF  (S( J+l) •GT.l  )  SI  J)  =  S( J  +  l)/2 
800  CONTINUE 
AL-BS 
J  J  =  0 

DO  900  A= 1 , AL 
DO  900  B  =  A i BL ,  B  S 
DO  900  C  =  B  ,  CL  » C  S 
DO  900  C-C  »  DL , DS 
DO  900  E  =  D»  EL ,  E  S 
DO  900  F  =  E  ,  F  L  ,  F  S 
DO  900  G-F , GL , GS 
DO  900  H=G»HL,HS 
DO  900  I=Hf IL, IS 
DO  900  J=I,JL,JS 
DO  900  K= J , KL , KS 
DO  900  L=K,LL,LS 
DO  900  M=LtMLfMS 
J  J  —  J  J  "3"  1 

IF  (JJ.LE.M)  GO  TO  900 
T=X { J J ) 

X(JJ)  =  XU  ) 

X ( M ) -T 
T=Y { J J ) 

Y  {  J  J  )  =  Y  {  ^  ) 

Y { M  )  =T 

900  CONTINUE 
RETURN 
END 


SUBROUTINE  RECOVR  I L0G2N ,  RE,  IM) 


RECOVER?  THE  REAL  AND  THE  IMAGINARY  PARTS  OF  THE  FOURIER 
TRANS  FO  F  M  WHICH  HAS  BEEN  STORED  IN  COMPACT  HERMET I  AN  FORM 


r>  o  O 
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C 

c 

c 

c 

c 

c 

c 

c 


c 


THE  COMPACT  STORED  FCIRIER  TRANSFORM  IS  ORIGINALLY  IN  RE  C  ) 
SEE  TWO  R  FT  FOR  DETAILS. 

IN  EFFECT  THIS  SUBROUTINE  DOUBLES  THE  STORAGE  AREA  REQUIRED 
BUT  GIVES  EASE  IN  THE  MANIPULATION  OF  THE  COEFFICIENTS  WHERE 
FT XI  J)  =  RE  t  J)  -J-  I M  <  J) 

REAL  RE  UO)f  IM  (10) 

INTEGER  LOG2N  p  N,  N  UPGN2»  MORE 

N  =  2**L CG2N 
N  UPON  2  =  N/2 
MORE  =  N>  UPON  2-8-2 
I  M ( 1 )  =  0«0 
DO  10  J~  2  p  N  UPON  2 
I  =  N-J+2 
I  M  {  J  )  =  R  E  {  I  ) 

10  CONTINUE 

J  =  N  UPON  2  4-1 
I M { J )  =  C.O 
DO  20  J— MOREjN 
I  =  N-J  4-2 
I M ( J )  =  -RE { J ) 

RE ( J )  r  RE ( I ) 

20  CONTINUE 

RETURN 
END 


SUBROUTINE  ONE  C  FT  (L0G2N,X) 
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ONE  COMPLEX  FOURIER  TRANSFORM 

INTEGER  L0G2N 
REAL  X  (ID) 

C 

INTEGER  J  ,  JN»K,KN? N, N  OVER  2 
REAL  AIj/5R»ARGpBI?BRtCjPI  t  S  ?  T 
C 

N=2**( L0C-2N-1 ) 

P I  =  3  «  1 4 1  -  92654 

N  OVER  2  =  N / 2  + 1 

DO  100  J=2fN  OVER  2 

K=N+2-J 

JN= J+N 

KN=K+N 

AR=X( J)+> (K) 

AI=X(KN)-X{JN) 


. 


' 


BR=X I J ) -X { K } 

Bl=X{ JN) +XCKN) 

ARG=PI*FLOAT( J-l 3/FL0ATIN ) 

C-CCS ( ARG ) 

S=SIN( ARG) 

T=BR*C+B  l*S 
BI=BI*C-BR*S 
BR  =  T 

X( J )  =  AR—  8  I 
X ( K }=AR+BI 
X( JN)=BR+AI 
X(KN)=BR-AI 
100  CONTINUE 

A  R  =  X  i  1 ) *X(N+1 ) 

BR=  X { 1 )  —  X 1 N  +  l ) 

X( 1 )=AR 
X(N  +  D— 8R 

CALL  MR  ID  FT  ( LOG 2N- 1 , X { 1 ) » X ( N  + 1 ) ) 
CALL  R  SUB  B  0  (  L0G2N-1 ,  X  <  1  )  ) 

CALL  R  SUB  B  0  ( L0G2N-1 , X ( N+l ) } 

CALL  R  SUB  B  0  <L0G2N,X) 

RETURN 

END 


